Implementing Linear Combination of Unitaries on Intermediate-term Quantum Computers

Shantanav Chakraborty
CQST and CSTAR, International Institute of Information Technology Hyderabad, Telangana, India

Over the years, the framework of Linear combination of unitaries (LCU) has been extremely useful for designing a plethora of quantum algorithms. In this work, we explore whether this widely applicable paradigm can be implemented on quantum computers that will be available immediately after the current NISQ stage. To this end, we develop three variants of LCU and apply each, to quantum algorithms of practical interest.

First, we develop a physically motivated, continuous-time analogue of LCU (“Analog LCU”). This technique, implementable on hybrid qubit-qumode systems, is simpler than its discrete-time counterpart. We use this method to develop analog quantum algorithms for ground state preparation and quantum linear systems. We also develop a randomized quantum algorithm to sample from functions of Hamiltonians applied to quantum states (“Single-Ancilla LCU”). This approach repeatedly samples from a short-depth quantum circuit and uses only a single ancilla qubit. We use this to estimate expectation values of observables in the ground states of a Hamiltonian, and in the solution of quantum linear systems. This method is suitable for early fault-tolerant quantum computers. Our third approach stems from the observation that for several applications, it suffices to replace LCU with randomized sampling of unitaries according to the distribution of the LCU coefficients (“Ancilla-free LCU”). This is particularly useful when one is interested in the projection of a quantum state implemented by an LCU procedure in some subspace of interest. We demonstrate that this technique applies to the spatial search problem and helps establish a relationship between discrete and continuous-time quantum walks with their classical counterparts. Our work demonstrates that generic quantum algorithmic paradigms, such as LCU, can potentially be implemented on intermediate-term quantum devices.

§ INTRODUCTION

We are currently in an era of quantum computing where theoretical advancements have been accompanied by drastic improvements in experimental capabilities [1, 2, 3, 4, 5]. With rapid progress being made, it is reasonable to envision a stage in the near future, where quantum computing will transition away from the NISQ era [6, 7]. Quantum devices, available immediately after the current NISQ stage, will most likely not have the capabilities of a large-scale, fully-programmable, fault-tolerant quantum computer. These devices would have short circuit depth and only a limited number of logical qubits, the so-called early fault-tolerant quantum computers [8, 9, 10, 11, 12]. On the other hand, for particular quantum technological platforms, it might be possible to engineer certain specific interactions more precisely, and for longer time-scales than others. For instance, it might be easier to engineer hybrid qubit-qumode systems in the intermediate-term  [13, 14, 15, 16, 3, 17, 18] as many of the most promising quantum technological platforms such as superconducting systems [13], ion-traps [19], and photonic systems [20], naturally have access to continuous variables. We refer to such devices, which will become available shortly after the current stage, as “intermediate-term quantum computers”.

It is thus crucial to develop quantum algorithms of practical interest that are implementable on intermediate-term quantum computers. Indeed, quantum algorithms tailored to early fault-tolerant quantum computers are already being developed [8, 9, 10, 11, 12]. With many quantum technological platforms vying for supremacy, it is  also essential to develop physically motivated quantum algorithms that can exploit the degrees of freedom that are naturally available to such platforms.

Over the years, only a few quantum algorithmic frameworks have found broad applicability: they can be used to solve various problems of interest. However, such frameworks are only implementable on fully-fault tolerant quantum computers, which might be decades away. In this work, we ask the following question:

Can we bring generic quantum algorithmic paradigms closer to implementation on intermediate-term quantum devices?

The framework of Linear Combination of Unitaries (LCU) is one such paradigm. Over the years, it has been widely applied and has been central to the development of a plethora of useful quantum algorithms ranging from Hamiltonian simulation [21, 22, 23, 24], quantum linear systems [25, 26] and differential equations [27, 28, 29], quantum walks [30, 31, 32], ground state preparation [33, 34, 35] and a large-class of optimization problems [36, 37]. In this work, we significantly enhance the applicability of the LCU framework to intermediate-term quantum devices: namely hybrid quantum systems and  early fault-tolerant quantum computers.

Given a Hermitian matrix $H$, LCU implements any function $f(H)$ that can be well-approximated by a linear combination of unitaries. Despite its broad applicability, LCU has its drawbacks when it comes to being implementable in the intermediate-term. First, for many problems of interest, there is a significant overhead in terms of the number of ancilla qubits needed. Additionally, the procedure requires implementing highly controlled-unitary operations, which can be challenging for near to intermediate-term quantum computers. In this work, we demonstrate that the issues concerning LCU can be addressed, and the framework can be used to develop new quantum algorithms that are more amenable to implementation. To this end, we introduce three new approaches to implementing LCU and demonstrate that each can be applied to develop quantum algorithms of practical interest.

Firstly, we develop “Analog LCU”, a physically motivated, continuous-time analogue for implementing a linear combination of unitaries. This technique requires coupling the system Hamiltonian $H$ to a continuous-variable ancilla system (such as a one-dimensional quantum Harmonic oscillator), initialized in some easy-to-prepare continuous-variable quantum state (such as a Gaussian). The overall system is then evolved according to the resulting interaction Hamiltonian. Although this approach requires a continuous-variable ancilla register, the overall algorithm is simpler. Moreover, this technique might be particularly useful for intermediate-term quantum computers as such interactions can already be engineered on several quantum technological platforms. Examples of discrete systems coupled to continuous-variable ones include ion traps and superconducting systems [13, 14, 16, 15, 17]. We show that this approach can be used to develop novel analog quantum algorithms for ground state preparation and solving quantum linear systems.

Our second approach is a randomized quantum algorithm that samples from functions of Hamiltonians applied to quantum states. More precisely, for any initial state $\rho_0$, and observable $O$, we estimate $\Tr[Of(H)\rho_0f(H)^{\dag}]$ to arbitrary accuracy. This technique, which we refer to as “Single-Ancilla LCU”, is a generalization of the Hamiltonian simulation algorithm of [9], and is implementable on early fault-tolerant quantum computers. The hallmark of quantum algorithms run on such devices is that they minimize the amount of quantum resources consumed [8, 12]. Typically, such quantum algorithms repeatedly sample from a short-depth quantum circuit, followed by classical post-processing. Our randomized quantum algorithm requires only one ancilla qubit that acts as a control and requires implementing two (controlled) unitaries sampled according to the distribution of the LCU coefficients, followed by a single-shot measurement. By repeating this simple quantum circuit, one obtains samples whose average converges to the expectation value we seek to estimate. Another advantage of this method is in the case where $H$ can itself be expressed as a linear combination of unitaries, i.e. $H=\sum_{j} c_j H_j$, where $H_j$ is unitary. Then in order to implement any $f(H)$ that can be well-approximated by some Fourier series, one can also take advantage of the recently developed randomized Trotter methods [38, 39, 9, 40]. The overall decomposition of $f(H)$ would be a linear combination of products of $H_j$, which can be repeatedly sampled (according to the distribution of the resulting LCU coefficients) while adding no overhead in terms of ancilla qubits. These features also make it appealing for early fault-tolerant quantum computers. We apply this method to estimate the expectation values of observables in the ground states of a Hamiltonian and the solution of quantum linear systems.

For several applications, it is sufficient to obtain only a projection of the LCU state $f(H)\ket{\psi_0}$ in some subspace. In such scenarios, we show that the ancilla registers can be dropped entirely. We call this the “Ancilla-free LCU” technique. This approach involves randomly sampling unitaries, according to the distribution of the LCU coefficients without any ancilla registers. More precisely, for any matrix $H$ if $f(H)\approx \sum_{j}c_j U_j$, we sample $U_j$ according to the distribution $\mathcal{D}\sim \{c_j/\norm{c}_1\}$. This results in some average density matrix $\rho$, for which the projection in this subspace can be proven to be at least as large. This is typically the case for quantum walk-based problems, where we are interested in obtaining a state that has a good overlap with some subset of nodes in a graph. For instance, consider quantum spatial search algorithms – where we are interested in preparing a state that has a good overlap with the marked nodes of any ergodic, reversible Markov chain. Consequently, we use this technique to design optimal spatial search algorithms by discrete-time and continuous-time quantum walks, also placing recent results in this context [31, 41, 32]. Note that for these algorithms, we do not need any extra registers (other than the walk space). In addition to this, these techniques allow us to connect discrete and continuous-time quantum walks with their classical counterparts. We believe that this approach can be applied to a wide class of problems, beyond quantum walks, such as quantum simulated annealing [42] and quantum metropolis algorithms [43, 44].

Ever since the work of Childs [45], one of the long-standing open problems has been to obtain a relationship between discrete-time and continuous-time quantum walks. Childs showed that one can obtain a discrete-time quantum walk for any Hamiltonian. However, given access to the evolution operator of a continuous-time quantum walk, it is not known whether one can obtain a discrete-time quantum walk. Virtually no progress has been made in this direction in decades. Along the way, using quantum singular value transformation (QSVT) [46, 47], we show how one can obtain discrete-time quantum walks from continuous-time quantum walks (and vice versa), thereby making significant progress on this problem.

The paper is organized as follows. We begin by providing a brief overview of the main results and also relate our work to prior results in Sec. <ref>. In Sec. <ref>, we review some basic definitions and techniques that we will be using in this article. We formally describe the three different approaches to implementing LCU in Sec. <ref>. The rest of the article involves applying these techniques to develop new quantum algorithms. In Sec. <ref>, we make use of our techniques to develop new quantum algorithms for ground state preparation of Hamiltonians and also to sample from their ground states. We develop analog quantum linear systems algorithms and also show how to sample from the solution of quantum linear systems in Sec. <ref>. The quantum algorithms developed in Sec. <ref> and Sec. <ref> mainly make use of the “Analog LCU” and the “Single-Ancilla LCU” frameworks, respectively. In Sec. <ref>, we make use of the “Ancilla-free LCU” technique to relate discrete and continuous-time quantum walks with their classical counterparts. In addition, we also show how one can obtain discrete-time quantum walks from continuous-time quantum walks and vice versa. Finally, we conclude and discuss open problems in Sec. <ref>.

§ SUMMARY OF OUR RESULTS

In this section, we outline the main results of this article. The LCU method involves implementing functions of matrices that can be approximated by linear combinations of unitaries. We begin by briefly outlining each of the three variants of implementing LCU that we develop in this article. Then we state the complexities of the various algorithms to which we apply these techniques. Finally, we also state some other related results that we develop along the way.

Summary of the main results – The three approaches to LCU and their applications

§.§ Analog LCU: Coupling discrete systems with continuous variable systems

Firstly, we develop a more physical model for LCU in continuous-time. Consider any Hamiltonian $H$, consider any $f(H)$ that can be well approximated by a truncated Fourier transform, i.e.,

\norm{f(H)-\int_{a}^{b}~dz~c(z) \cdot e^{-iHzt}}\leq \varepsilon,

where $c:\mathbb{R}\mapsto \mathbb{R}\backslash \{0\}$. Then by a purely continuous-time procedure, for any initial state $\ket{\psi_0}$, we can prepare a state that $O\left(\varepsilon/\norm{c}_1\right)$-close to $\frac{f(H)\ket{\psi_0}}{\norm{f(H)\ket{\psi_0}}}$,
where $\norm{c}_1=\int_{a}^{b} dz~c(z)$. This approach requires coupling $H$ to a continuous variable ancilla system (such as a one-dimensional quantum Harmonic oscillator) prepared in a continuous variable state. We show in this work that for several applications, this state is easy to prepare (such as a Gaussian). The overall system is then evolved according to the interaction Hamiltonian $H'=H\otimes \hat{z}$ for an appropriate time $T$. This approach is not only simpler than its discrete-time counterpart, but such hybrid qubit-qumode coupling can be implemented in a number of quantum technological platforms such as trapped ions, Cavity (or Circuit QED), photonic systems, and superconducting qubits [13, 14, 16, 15, 17]. Our motivation is not only to provide an alternative approach to implementing LCU but make this paradigm more implementable for intermediate-term quantum devices.

Applications:  We apply these techniques to develop analog quantum algorithms for ground state preparation and for solving quantum linear systems. As mentioned previously, our aim is to couple the system Hamiltonian with an ancillary continuous-variable system.

(a) Ground state preparation:  Given a Hamiltonian $H$, we couple this system with a continuous variable ancillary system via the interaction Hamiltonian $H'=H\otimes \hat{z}$. The ancilla system is prepared in an easy-to-prepare continuous variable state, namely a Gaussian. We show that given an initial state $\ket{\psi_0}$ that has an overlap of at least $\eta$ with the ground state, simply evolving the system according to $H'$ results in a state proportional to $f(H)\ket{\psi_0}$ in the first register, where $f(H)=e^{-tH^2}$. We show that, with probability $\eta^2$, this state is $\varepsilon$-close to the ground state of $H$ (provided its ground energy is known up to some precision). The overall time required is

where $\Delta$ is the spectral gap of $H$ (Lemma <ref>). This quantum algorithm appeared in [32] and also independently in Ref. [35]. Here we place this in the context of “Analog LCU” it provides useful intuition for (i) the quantum linear systems algorithms we develop using similar techniques and (ii) the “Single-Ancilla LCU”method to solve this problem.

(b) Quantum linear systems:  We provide two quantum algorithms for this problem. For both these problems, we couple $H$ to two ancillary continuous variable systems (Harmonic oscillators), i.e. $H'=H\otimes \hat{y}\otimes \hat{z}$.  The first approach works for any Hermitian matrix $H$ with eigenvalues in the domain $[-1,-1/\kappa]\cup [1/\kappa,1]$, where $\kappa$ is an upper bound on the condition number (ratio between the maximum and the minimum non-zero eigenvalue) of $H$. The first register contains the initial state $\ket{b}$, the second register is prepared in the first excited state of the first Harmonic oscillator, while the third register is prepared in the ground state of a “particle in a ring” of unit radius [48].

This algorithm (see Sec. <ref>) can be seen as an analog variant of the quantum linear systems algorithm of Childs, Kothari and Somma [25]. In order to obtain a quantum state that is $\varepsilon/\kappa$-close to $\ket{x}=\frac{H^{-1}{\ket{b}}}{\norm{H^{-1}{\ket{b}}}}$ in the first register, with overlap at least $1/\kappa$, we require evolving the system according to $H'$ for a time

Typically in continuous variable systems, Gaussian states are easier to prepare, and engineer [49]. Thus, we also provide an analog quantum algorithm for solving quantum linear systems (for positive semidefinite Hamiltonians) in which both the ancilla registers are now prepared in Gaussian states. Evolving this system according to $H'$ prepares a state that is $(\varepsilon/\kappa)^{3/2}$ - close to $\ket{x}$, with overlap $\Omega(1/T)$ in time

Although the complexity is worse than the first approach, this quantum algorithm requires preparing only Gaussian states, which we expect to be easier for intermediate-term quantum computers to implement.

§.§ Single-Ancilla LCU: Controlled application of randomly sampled unitaries followed by single-shot measurement

Given any Hamiltonian $H$, we develop a randomized quantum algorithm that estimates expectation values $\Tr[Of(H)\rho_0 f(H)^{\dag}]$, to arbitrary accuracy, where $O$ is any observable and $f(H)$ can be well-approximated by a linear combination of unitaries. The quantum algorithm is suitable for implementation on early fault-tolerant quantum computers: it repeatedly runs the simple short-depth quantum circuit shown in Fig. <ref> and requires only one ancilla qubit.

$\ket{+}$        1     1     [2]X⊗O

$\rho_0$            V_1     V_2

Quantum circuit corresponding for the “Single-Ancilla LCU” procedure. For $f(H)\approx \sum_{j}c_j U_j$, repeated runs of this short-depth quantum circuit can estimate $\Tr[O f(H)\rho_0 f(H)^{\dag}]$, to arbitrary accuracy. For this, $V_1$ and $V_2$ are sampled at random according to $\mathcal{D}\sim \{c_j/\norm{c}_1, U_j\}$. Each run of the circuit outputs a random variable corresponding to the outcome of the measurement of the observable $X\otimes O$. Overall we need to repeat this circuit $T$ times, which is large enough for the sample mean of the $T$ observations to converge to the expectation value.

Our approach is a generalization of the Hamiltonian simulation procedure of Faerhmann et al. [9], wherein the authors used this circuit to generate randomized multi-product formulas. We generalize this to implement any $f(H)$ which can be approximated by an LCU. The quantum circuit applies controlled and anti-controlled versions of $V_1$ and $V_2$, respectively where $V_1$ and $V_2$ are sampled according to $ \{c_j/\norm{c}_1, U_j\}$. This is followed by a measurement of the observable $X\otimes O$ which gives the output of a single run.  This procedure is then repeated enough times so that the sample mean of the outcomes converges to $\Tr[Of(H)\rho_0 f(H)^{\dag}]$. The overall procedure is formally stated in Algorithm <ref> whose correctness we prove via the following theorem:
[Sampling from functions of Hamiltonians applied to quantum states]theoremthmbodyrandomizedtimeevolve

Let $\varepsilon, \delta, \gamma \in (0,1)$ be some parameters. Let $O$ be some observable and $\ket{\psi_0}$ be some initial state. Suppose there is a Hermitian matrix $H\in\mathbb{R}^{N\times N}$, such that $\norm{f(H)-\sum_{j}c_j U_j}\leq \gamma$, where $U_j$ is some unitary such that

\gamma \leq \dfrac{\varepsilon}{6\norm{O}\norm{f(H)}}.

Furthermore, let

T\geq \dfrac{8\norm{O}^2\ln(2/\delta)\norm{c}_1^4}{\varepsilon^2}.

Then, Algorithm <ref> estimates $\mu$ such that

\left|\mu - \Tr[Of(H)\ket{\psi_0}\bra{\psi_0}f(H)^{\dag}]\right| \leq \varepsilon,

with probability at least $1-\delta$, using of one ancilla qubit and $T$ runs of the quantum circuit in Fig. <ref>.

Overall, $T$ samples are required, and if $\mathcal{C}_j$ is the cost implementing $U_j$ and the state $\ket{\psi_0}$ can be prepared in cost $\mathcal{C}_{\psi_0}$, then the cost of each run of the quantum circuit is upper bounded by $\max_j C_j + \mathcal{C}_{\psi_0}$. Furthermore, when $f(H)$ has a well-defined Fourier series, and each $U_j=e^{-ijH}$,  the cost of each run is determined by the cost of the underlying Hamiltonian simulation algorithm. For instance, consider that $H=\sum_{j} \alpha_j \Lambda_j$, where $\Lambda_j$'s are local unitary operators (a string of Pauli matrices). Then one can use any of the recently developed randomized Trotter-based approaches for Hamiltonian simulation such as qDRIFT [38, 39] or randomized multi-product formulas [50, 9], and others [51, 40]. These methods do not require ancilla qubits and could be easily incorporated into our algorithm.

Applications:  We use this technique to sample from the ground states of Hamiltonians and also from the solution of quantum linear systems.

(a) Ground state preparation:  Suppose $H$ is a  Hamiltonian with ground state $\ket{v_0}$ and spectral gap $\Delta$, with so that its ground energy is known to a certain precision. Also, consider some initial state $\ket{\psi_0}$ that can be prepared with cost $\tau_{\psi_0}$ and has an overlap of at least $\eta$ with $\ket{v_0}$. For any observable $O$, by considering the LCU of the function $f(H)=e^{-tH^2}$, we can use Algorithm <ref> to output $\mu$ such that

\left|\mu - \braket{v_0|O|v_0}\right| \leq \varepsilon,

with probability at least $1-\delta$, using

T\geq O\left(\dfrac{\norm{O}^2\ln(2/\delta)}{\varepsilon^2\eta^4}\right),

runs of a quantum circuit where the cost of each run is $\tau_{\max}+\tau_{\psi_0}$, where

\tau_{\max}=O\left(\dfrac{1}{\Delta}\log\left(\dfrac{\norm{O}}{\varepsilon\eta}\right)\right).

This has been formally stated in Theorem <ref>.

(b)  Quantum linear systems:  Suppose we have a Hermitian matrix $H$ such that its eigenvalues in $[-1, -1/\kappa]\cup [1/\kappa, 1]$. Let us assume that the initial quantum state $\ket{b}$ can be prepared in cost $\tau_b$. Then, we show that we can use Algorithm <ref> to output some

\left|\mu-\braket{x|O|x}\right|\leq \varepsilon,

with probability at least $1-\delta$, using

T\geq O\left(\dfrac{\norm{O}^2\kappa^4\log^2\left(\dfrac{\norm{O}\kappa}{\varepsilon}\right)\ln(2/\delta)}{\varepsilon^2}\right),

runs of the quantum circuit in Fig. <ref> where the cost of each run is $\tau_{\max}+\tau_b$, where

\tau_{\max}=O\left(\kappa\log\left(\dfrac{\norm{O}\kappa}{\varepsilon}\right)\right).

This approach makes use of the LCU decomposition of $f(H)=H^{-1}$ in Ref. [25] (See Theorem <ref>).
We also develop an alternative approach that has a slightly better sample complexity (without the log factor). This makes use of the key insight of the recently developed adiabatic approaches for solving quantum linear systems [52, 53, 54, 55]. Namely, we can construct some Hamiltonian $H'$, given any $H$ such that the $0$-eigenstate of $H'$ corresponds to the solution of the quantum linear systems, i.e. $\ket{x}=H^{-1}\ket{b}/\norm{H^{-1}\ket{b}}$. So, we are able to exploit this connection to use the procedure in (a) above (For details see Theorem <ref>). However the construction of $H'$ requires an additional qubit.

§.§ Ancilla free LCU: Randomized sampling of unitaries

Suppose for some Hermitian matrix $H$, we intend to implement $f(H)$ such that

\norm{f(H)-\sum_{j=1}^{M}c_j U_j}\leq \gamma,

where $\gamma\in [0,1)$, $c_j\in \mathbb{R}$ and $U_j$ is some unitary. For instance, it may be the case that $U_j=e^{-ijH}$.  We show that dropping the ancilla register altogether and simply sampling $U_j$ to $\mathcal{D}\sim \{c_j/\norm{c}_1\}$, allows us to prepare the average density matrix

\rho = \dfrac{1}{\norm{c}^2_1} \sum_{j=1}^{M} c_j U_j \ket{\psi_0}
\bra{\psi_0} U^{\dag}_j,

where $\norm{c}_1=\sum_{j=1}^{M} c_j$. In Sec. <ref>, we show that this suffices, for instance, if we are interested in the projection of $f(H)\ket{\psi_0}$ in some subspace of interest. We may instead prepare $\rho$, by sampling $U_j$ according to $\mathcal{D}$ (or evolving the initial state under $H$ for a random time). The projection of $\rho$ in this subspace is guaranteed to be at least as large. Formally, we prove the following theorem

[Randomized unitary sampling]theoremthmbodyunisamp

Let $H\in \mathbb{R}^{N\times N}$ is a Hermitian matrix. Also let $\varepsilon\in (0,1)$ and suppose $f:\mathbb{R}\mapsto\mathbb{R}$ be some function such that

\norm{f(H)-\sum_{j=1}^{M} c_j U_j}\leq \dfrac{\varepsilon}{3\norm{f(H)}},

for some unitaries $U_j$ and $c_j\in \mathbb{R}\backslash\{0\}$. For some initial state $\rho_0$, define the density matrix

\rho = \dfrac{1}{\norm{c}^2_1} \sum_{j=1}^{M} c_j U_j \rho_0 U^{\dag}_j,

where $\norm{c}_1=\sum_{j=1}^{M} c_j$.  Then, for any projector $\Pi$,

\norm{c}^2_1 \Tr\left[\Pi\rho\right]\geq \Tr[\Pi f(H)\rho_0 f(H)^{\dag}] -\varepsilon.

We prove Theorem <ref> in Sec. <ref>. Let us, for simplicity, let us assume that $\norm{c}^2_1 < 1$. Then the result of Theorem <ref> can be interpreted as follows: if we are interested in the projection of $f(H)\ket{\psi_0}$ in some subspace of interest, instead of implementing the LCU procedure, we can prepare $\rho$ by random unitary sampling. The projection of $\rho$ in this subspace is guaranteed to be is at least as large. The same approach would also work in the continuous-time setting described previously. This makes this technique suitable for application to the spatial search problem, where we are interested in preparing a quantum state with a good overlap with the space spanned by the marked vertices. It can also be applied to quantum simulated annealing [42], where the goal is to prepare a state that has a good overlap with the ground state. In this work, we apply this procedure to the former problem.

Given any ergodic, reversible Markov chain $P$ of $|X|=n$ nodes (out of which $M$ are marked), the spatial search algorithm finds a marked vertex. We discuss two separate quantum algorithms by discrete-time quantum walks that solve this problem quadratically faster than classical random walks. The first relies on fast-forwarding discrete-time  random walks and formalizes an unproven observation in [41]. The second quantum algorithm fast-forwards continuous-time random walks. Furthermore, we also briefly discuss the optimal spatial search algorithm by continuous-time quantum walk of Ref. [32]. The running time of all these algorithms scales as the square root of the hitting time of classical random walks (up to log factors), which is optimal. We demonstrate that the “Ancilla-free LCU” method captures all these quantum algorithms.

We briefly describe the results we obtain in both the discrete and the continuous-time settings:

(a) Discrete-time quantum walks:  We use the fact that if any Hamiltonian $H$ is encoded in the top-left block of a unitary $U_H$, we can obtain a block encoding [56, 26] of $H^t$ or $e^{-t(I-H)}$ by implementing a linear combination of Chebyshev polynomials of $H$. Both these procedures require roughly $O(\sqrt{t})$ cost to be implemented. When $H=D$ (the discriminant matrix of $P$), $D^t$ results in a fast-forwarding of discrete-time random walks [30]. One the other hand, when $H=e^{-t(I-D)}$, this is a fast-forwarding of continuous-time random walks. Now for the spatial search problem, we do not need to implement a full LCU procedure and can make use of Theorem <ref> instead.

Using the framework of interpolated Markov chains (see Sec. <ref> for details of these terms), and for a specific initial state $\ket{\sqrt{\pi_U}}$ (related to the stationary distribution of the interpolated random walk), we can show for the first spatial search algorithm, the “Ancilla-free LCU” procedure prepares an average density such that

\Tr[(I\otimes\Pi_M)\rho]\geq  \norm{\Pi_M D(s)^T\ket{\sqrt{\pi_U}}}^2-\varepsilon,

where $\Pi_M$ is a projection on to the marked subspace. In Ref. [31], the authors proved that the RHS of the aforementioned inequality is $\tilde{\Omega}(1)$, for $T=\widetilde{O}\left(\sqrt{HT}\right)$, for some randomly chosen value of $s\in [0,1)$ and $HT$ is the hitting time of a random walk on $P$ (See Algorithm <ref>). This also formalizes the observation of Ref. [41].

For our second spatial search algorithm by discrete-time quantum walk (See Algorithm <ref>), Theorem <ref>, prepares an average density matrix $\rho$ such that

\Tr[(I\otimes \Pi_M)\rho]\geq \norm{\Pi_M e^{T(D(s)-I)}\ket{\sqrt{\pi_U}}}^2-\varepsilon,

where again, the RHS can be proven to be in $\widetilde{\Omega}(1)$, for $T=\widetilde{O}\left(\sqrt{HT}\right)$.

(b) Continuous-time quantum walks:  The optimal quantum spatial search algorithm by continuous-time quantum walk [32] is yet another demonstration of ancilla-free LCU, where a quadratic speedup for this problem can be obtained simply by randomized time-evolution. For this, the key idea is that there exists a Hamiltonian $H$ that corresponds to a quantum walk on the edges of the underlying Markov chain $P$ such that $H^2\ket{\psi}\ket{\bar{0}}=[(D^2-I)\otimes I]\ket{\psi}\ket{\bar{0}}$. Then by using the analog LCU procedure to implement $f(H)=e^{-tH^2}$, one can implement $e^{t(D^2-I)}$. As before, for the spatial search problem we can bypass the LCU itself and use Theorem <ref> instead to obtain $\rho$ such that,

\Tr[(I\otimes \Pi_M)\rho]\geq \norm{\Pi_M e^{T(D^2(s)-I)}\ket{\sqrt{\pi_U}}}^2-\varepsilon,

where in [32], the RHS is $\widetilde{\Omega}(1)$, for $T=\widetilde{O}\left(\sqrt{HT}\right)$.

These algorithms allow us to relate discrete-time random walks and continuous-time random walks with their quantum counterparts.

§.§ Other results on quantum walks

Other than the aforementioned contributions, we also develop several additional results along the way. A long-standing open problem in this area has been to establish a relationship between discrete-time quantum walks and continuous-time quantum walks. In a seminal work, Childs showed that any dynamics generated by a continuous-time quantum walk can be simulated by a discrete-time quantum walk [45]. However, whether discrete-time quantum walks can be obtained from the continuous-time quantum walk evolution operator, has been unknown. In this work, we also make progress towards answering this question. In fact, we complete the picture by also connecting random and quantum walks in both discrete and continuous-time settings.

Summary of other results – Relationship between discrete and continuous-time quantum walks, and their classical counterparts.

We first show that given a block encoding of any Hamiltonian $H$, it is possible to obtain a discrete-time quantum walk on $H$. This was observed also in [41]. As we have already discussed, from this block encoding, it is possible to fast-forward discrete-time and continuous-time random walks. Also, following the result of [32], we know that a continuous-time quantum walk can fast-forward a continuous-time random walk.

* From discrete-time quantum walks to continuous-time quantum walks:  In this work, we also show that from any given unitary $U_H$, which is a block encoding of $H$ implemented in cost $T_H$, one can construct a block encoding of the Hamiltonian $H_P=i[U,\Pi_0]$ in cost $O(T_H)$, where $\Pi_0=I\otimes\ket{\bar{0}}\bra{\bar{0}}$ such that $\ket{\bar{0}}$ is some reference state. When $H$ corresponds to the discriminant matrix of an ergodic, reversible Markov chain $P$, $e^{-iH_Pt}$, corresponds to a continuous-time quantum walk on the edges of $P$. Thus, using Hamiltonian simulation algorithms, we are able to obtain a continuous-time quantum walk, starting from a discrete-time quantum walk.

* From continuous-time quantum walks to discrete-time quantum walks:  Suppose $H$ is some Hamiltonian that encodes the connectivity of some graph $G$. Then $U=e^{iH/2}$ is a continuous-time quantum walk on $G$. Using quantum singular value transformation, we explicitly show that given access to $U$, we can construct a $(1,3,\varepsilon)$ block encoding of $H$ in cost $O\left(\frac{1}{\varepsilon}\log(1/\varepsilon)\right).$ This allows us to obtain a discrete-time quantum walk from a continuous-time quantum walk. We discuss the subtleties of this approach and also possible improvements. The overall relationship between the different quantum and random walk frameworks is indicated in Fig. <ref>.

§.§ Prior work

In this section, we briefly sketch relevant prior work and relate them to the results we obtain. The linear combination of unitaries technique was first developed by Childs and Wiebe [21] to develop quantum Hamiltonian simulation algorithms based on multi-product formulas. Since then, LCU has been extensively used to develop improved quantum algorithms for Hamiltonian simulation [22, 24, 23]. Subsequently, LCU has been used to develop a wide variety of quantum algorithms for linear algebra, such as for solving quantum linear systems, and linear regression [25, 26], preparing ground states of Hamiltonians [33] and solving optimization problems [36, 37]. Many of these quantum algorithms are unified by the framework of quantum singular value transformation (QSVT) [46], which implements polynomial transformations to the singular values of a matrix. Although QSVT provides near-optimal query complexities for these problems and requires fewer ancilla qubits than LCU, the framework itself is not likely to be implemented on early fault-tolerant quantum computers. In order to implement polynomial transformations to the underlying matrix, QSVT requires a sequence of controlled operations and the overall quantum circuit can be of large depth. It is also unclear whether the QSVT-based techniques can be modified to make it more amenable for near/intermediate-term implementation. This remains an important open problem.

The main contribution of this article is to demonstrate that the framework of implementing LCU can be modified so that this framework, which is also applicable to a wide variety of problems, is implementable on intermediate-term quantum computers. As discussed in the previous sections, we introduce three main variants of LCU.

The first technique is a continuous-time variant of LCU, which is more physical. The key idea is to couple discrete systems with continuous-variable systems. Such interactions have been explored in the context of quantum phase estimation, where the system Hamiltonian is coupled with a one-dimensional free particle, acting as the pointer variable - which is the so-called von Neumann measurement model [57, 58]. In Ref. [32], the continuous-time quantum walk Hamiltonian $H$ was coupled to a one-dimensional quantum Harmonic oscillator to implement $e^{-tH^2}$. This was a key ingredient of their spatial search algorithm by continuous-time quantum walk. In this work, we formalize this technique and show that it is more widely applicable and in fact, can serve as a continuous-time variant to any LCU-based quantum algorithm. We develop an analog variant of the quantum linear systems algorithm of Childs et al. [25] and a new quantum algorithm for this problem (using only Gaussian states) that is more suited for intermediate-term implementation.

The “Single-Ancilla LCU" technique is a generalization of the work by Faerhmann et al. [9], where
similar techniques were used for Hamiltonian simulation. We apply this technique to develop randomized quantum algorithms for ground state preparation and solving quantum linear systems that are suitable for implementation on early fault-tolerant quantum computers. For both these problems, several quantum algorithms have been developed over the years. The first quantum algorithm for ground state preparation involved using Hamiltonian simulation along with quantum phase estimation [59]. Subsequently, Refs. [33, 34] took advantage of the fact that functions of Hamiltonians can be expressed as a linear combination of unitaries to develop fast quantum algorithms for ground state preparation and estimation. A QSVT-based quantum algorithm has also been developed recently [60]. The problem of preparing ground states of Hamiltonians is considered to be one of the first problems to be solved on early fault-tolerant quantum computers. In this context, several quantum algorithms have been proposed for ground state preparation as well as for estimating properties of ground states [61, 12, 11].

Ever since the seminal algorithm by Harrow, Hassidim and Lloyd [62], quantum linear systems has been analyzed extensively. In particular, the LCU-based approach of Childs, Kothari, and Somma [25] provided a linear dependence on the condition number of the underlying sparse matrix and an exponentially improved dependence on the error. This algorithm was improved to also work for non-sparse matrices [63] and in the more general block encoding framework [64]. Recently, QSVT-based approaches for this problem have also been developed [46, 47, 64]. Another direction of research has been to develop quantum algorithms for this problem in the adiabatic quantum computing framework [52, 53, 55]. The possibility of applying quantum linear systems algorithm on the near-term quantum devices has been explored in [65].

Finally, we apply the “Ancilla-free LCU” technique to develop optimal quantum spatial search algorithms. Recently, LCU techniques have been used to develop several quantum walk algorithms. The quantum fast-forwarding scheme by Apers and Sarlette [30] quadratically fast-forwards the dynamics of a discrete-time random walk by implementing a linear combination of discrete-time quantum walk steps. Recently,  Ambainis et al. [31] proved that for the spatial search problem, fast-forwarding $T$ random walk steps on an interpolated Markov chain, prepares a quantum state that has a $\tilde{\Omega}(1)$ overlap with the marked space for $T=\widetilde{O}(\sqrt{HT})$, where $HT$ is the classical hitting time of the random walk. Thus, their LCU-based quantum spatial search algorithm for discrete-time quantum walks completely solves the spatial search problem quadratically faster than classical random walks, for any number of marked nodes. This closed a long line of work which made partial progress towards solving this problem. Subsequently, Apers et al. [41] provided a unified framework that connected the different variants of discrete-time quantum walk search. Therein, the authors observed that the LCU procedure of [31] could be replaced by randomly sampling quantum walk steps but no proof was provided. In the continuous-time quantum walk framework, whether the spatial search problem offered a generic quadratic speedup was also open for a long time, and was only recently solved in [32]. Their analog quantum algorithm indeed managed to bypass the LCU procedure by evolving the system under the quantum walk Hamiltonian for a random time.

Establishing a relationship between discrete-time quantum walks and continuous-time quantum walks from both directions has been a long-standing open problem in the field. Ever since Childs' seminal work in this direction [45], virtually no progress has been made. This also explains why efforts to prove the optimality of quantum spatial search were undertaken in parallel in the discrete-time and the continuous-time settings. We show how one can obtain discrete-time quantum walks from continuous-time quantum walks and vice versa, making significant progress on this problem.

In the following section, we briefly define some of the key concepts that we use to derive our results.

§ PRELIMINARIES

In this section, we introduce some of the topics that we deal with in this article as well as discuss the key algorithmic primitives required to develop our results. We begin by introducing the LCU framework.

§.§ Linear Combination of Unitaries

We will begin by stating the general framework of Linear Combination of Unitaries (LCU). Let us define an operator $V=\sum_{j=1}^{M}c_j U_j$,  such that each $U_j$ is a unitary matrix. Without loss of generality let us define the parameters $c_j\in \mathbb{R}\backslash \{0\}$. Note that these parameters can take complex values as well, but in such cases, we can always absorb its imaginary part into the definition of the corresponding $U_j$ itself. Note that $V$ is not necessarily unitary and the LCU approach allows us to implement any such $V$.

Let $R$ be a state preparation unitary such that

where $c=(c_0, c_1,\dots c_M)^T$. Furthermore, define the controlled unitary

W=\sum_{j=1}^{M} \ket{j}\bra{j}\otimes U_j.

Then, it is easy to verify that the state

\begin{equation}
\ket{\psi_t}=(R^{\dag}\otimes I) W (R\otimes I)\ket{\bar{0}}\ket{\psi}=\dfrac{1}{\norm{c}_1}\ket{\bar{0}}\sum_{j=1}^{M} c_jU_j\ket{\psi}+\ket{\Phi}^{\perp},
\end{equation}

where $(\ket{\bar{0}}\bra{\bar{0}}\otimes I)\ket{\Phi}^{\perp}=0$. Thus, by post-selecting of obtaining $\ket{\bar{0}}$ in the first register, we obtain the state $V\ket{\psi}/\norm{V\ket{\psi}}$ in the second register, with probability $\norm{V\ket{\psi}}^2/\norm{c}^2_1$. The dominant factor in the complexity of implementing this procedure is $C_j$ the cost of implementing controlled $U_j$. This is then upper bounded by $C_{\max}=\max_{j} C_j$.

Let us now explore the applicability of this procedure. Given any Hermitian matrix $H\in \mathbb{R}^{N\times N}$ with spectral decomposition $H=\sum_{j=1}^{N}\lambda_j\ket{v_j}\bra{v_j}$, define $f(H)=\sum_{j=1}^{N}f(\lambda_j)\ket{v_j}\bra{v_j}$. Now suppose $f(H)$ can be well approximated by linear combinations of unitaries. More precisely, for some $\varepsilon\in (0,1)$ suppose,

\norm{f(H)-\sum_{j=1}^M c_j U_j}\leq \varepsilon.

For instance, $f(H)$ could be well approximated by a Fourier series, in which case, $U_j=e^{-ijH}$. Since, this is indeed the case, for several functions $f(x)$, LCU provides a versatile framework to implement these matrix functions. Consequently, several near-optimal quantum algorithms have been designed in this framework that has wide applicability such as quantum algorithms for linear systems [25], ground state preparation [33, 34], sampling from thermal states [36, 66], Hamiltonian simulation [21, 22, 23, 24] and others. However, as observed from the construction of $B$ and $W$, these quantum algorithms will be implementable only on a full-scale fault-tolerant quantum computer.

§.§ Block encoding and Quantum Singular Value Transformation

For many of the algorithms in this work, we will consider the framework of block encoding, wherein it is assumed that the input matrix $H$ (up to some sub-normalization) is stored in the left block of some unitary. The advantage of the block encoding framework, which was introduced in a series of works [67, 26, 46], is that it can be applied to a wide variety of input models.

Suppose that $H$ is an $s$-qubit operator, $\alpha, \varepsilon \in \mathbb{R}^+$ and $a \in \mathbb{N} $, then we say that the $(s + a)$-qubit unitary $U_H$ is an $(\alpha, a, \varepsilon)$-block encoding of $H$, if

\begin{equation}
\norm{ H - \alpha (\bra{0}^{\otimes a} \otimes I)U_H(\ket{0}^{\otimes a} \otimes I) } \leq \varepsilon .
\end{equation}

Let $\ket{\psi}$ be an $s$-qubit quantum state. Then applying $U_H$ to $\ket{\psi}\ket{0}^{\otimes a}$ outputs a quantum state that is $\frac{\varepsilon}{\alpha}$-close to

\dfrac{H}{\alpha}\ket{\psi}\ket{0}^{\otimes a}+\ket{\Phi^{\perp}},

where $\left(I_s\otimes \ket{0}^{\otimes a}\bra{0}^{\otimes a}\right)\ket{\Phi^{\perp}}=0$. Equivalently, suppose $\tilde{H} := \alpha\paren{\bra{0}^{\otimes a} \otimes I_s} U_H \paren{\ket{0}^{\otimes a} \otimes I_s}$ denotes the actual matrix that is block-encoded into $U_H$, then $\norm{H-\tilde{H}}\leq\varepsilon$.

Quantum Singular Value Transformation (QSVT) applies a polynomial transformation to the singular values of a block-encoded matrix [46]. Formally, let $P \in \mathbb{C}[x]$ be a polynomial of degree $d \ge 2$, such that

* $P$ has parity-$(d\mod 2)$,

* $\forall x \in [-1, 1] : \abs{P(x)} \leq 1$,

* $\forall x \in (-\infty, -1] \cup [1, \infty) : \abs{P(x)} \geq 1$,

* if $d$ is even, then $\forall x \in \rr : P(ix) P^*(ix) \geq 1$.

Then, QSVT allows us to implement any polynomial $P(x)$ that satisfies the aforementioned requirements. Next, we formally introduce QSVT formally via the following theorem.

Suppose $A\in \mathbb{R}^{N\times d}$ is a matrix with singular value decomposition $A=\sum_{j=1}^{d_{\min}}\sigma_j\ket{v_j}\bra{w_j}$, where $d_{\min}=\min\{N,d\}$ and $\ket{v_j}$ $(\ket{w_j})$ is the left (right) singular vector with singular value $\sigma_j$. Furthermore, let $U_A$ be a unitary such that $A=\widetilde{\Pi}U_A\Pi$, where $\Pi$ and $\widetilde{\Pi}$ are orthogonal projectors. Then, for any QSP polynomial $P(x)$ of degree $n$, there exists a vector $\Phi=\left(\phi_1,\phi_2,\cdots\phi_n\right)\in\mathbb{R}^n$ and a unitary

\begin{equation}
\label{eqn:qsvt_sequence}
U_{\Phi} = \begin{cases}
e^{i \phi_1 ( 2 \widetilde{\Pi} - I )} U_A \left[ \prod_{k = 1}^{(n - 1) / 2} e^{i \phi_{2k} ( 2 \widetilde{\Pi} - I )} U_A^{\dagger} e^{i \phi_{2k+1} ( 2 \widetilde{\Pi} - I )} U_A \right], & n \text{ is odd } \\
\left[ \prod_{k = 1}^{n/ 2} e^{i \phi_{2k-1} ( 2 \widetilde{\Pi} - I )} U_A^{\dagger} e^{i \phi_{2k} ( 2 \widetilde{\Pi} - I )} U_A \right], & n \text{ is even},
\end{cases}
\end{equation}

such that

\begin{equation}
\label{eqn:qsvt_sequence_projected}
P^{SV} (A) = \begin{cases}
\widetilde{\Pi} U_{\Phi} \Pi, & n \text{ is odd} \\
\Pi U_{\Phi} \Pi, & n \text{ is even},
\end{cases}
\end{equation}

where $ P^{SV} (A) $ is the polynomial transformation of the matrix $ A $ defined as

\begin{equation}
\label{eqn:poly_sv_transformation_of_matrix}
P^{SV} (A) := \begin{cases}
\sum_j P(\sigma_j) \ket{v_j} \bra{w_j}, & P \text{ is odd} \\
\sum_j P(\sigma_j) \ket{w_j} \bra{w_j}, & P \text{ is even}.
\end{cases}
\end{equation}

Theorem <ref> tells us that for a $P$ of degree $n$, we can implement $P^{SV}(A)$ using one ancilla qubit, $\Theta(n)$ applications of $U_A$, $U^{\dag}_A$ and controlled reflections $I-2\Pi$ and $I-2\widetilde{\Pi}$. Furthermore, if in some well-defined interval, some function $f(x)$ is well approximated by an $n$-degree polynomial $P(x)$, then Theorem <ref> also allows us to implement a transformation that approximates $f(A)$, where

\begin{equation}
\label{eqn:sv_transformation_of_matrix}
f(A) := \begin{cases}
\sum_j f(\sigma_j) \ket{v_j} \bra{w_j}, & P\text{ is odd} \\
\sum_j f(\sigma_j) \ket{w_j} \bra{w_j}, & P \text{ is even}.
\end{cases}
\end{equation}

The following theorem from Ref. [46] deals with the robustness of the QSVT procedure, i.e. how errors propagate in QSVT. In particular, for two matrices $A$ and $\tilde{A}$, it shows how close their polynomial transformations ($P^{SV}(A)$ and $P^{SV}(\widetilde{A})$, respectively) are, as a function of the distance between $A$ and $\tilde{A}$.

Let $P \in \mathbb{C}[x]$ be a $d$-degree polynomial that satisfies the requirements of QSVT. Let $A, \tilde{A} \in \mathbb{C}^{N \times M}$ be matrices of spectral norm at most 1. Then,

\begin{equation*}
\norm{P^{SV}(A) - P^{SV}(\tilde{A})} \leq 4d \sqrt{\norm{A - \tilde{A}}}.
\end{equation*}

Having discussed the preliminary concepts, in the next section, we explain the three variants of LCU we consider in this article.

§ NEW APPROACHES FOR IMPLEMENTING LINEAR COMBINATIONS OF UNITARIES

In this section, we present the key technical contributions of this work. We begin by describing a purely continuous-time variant of the LCU technique.

§.§ Analog LCU: coupling a discrete primary system to a continuous-variable ancilla

Suppose we have some Hermitian matrix $H\in \mathbb{R}^{N\times N}$ of unit spectral norm and we wish to implement $f(H)$ for some function $f:[-1,1]\mapsto \mathbb{R}$ which satisfies:

\left|f(x)-\int_{a}^{b} dz~c(z)\cdot e^{-itxz}\right|\leq \varepsilon,

where $c:\mathbb{R}\mapsto\mathbb{R}^{+}\backslash\{0\}$. For instance, this can be a Fourier transform of $f(x)$, in which case $a=-\infty$, $b=+\infty$ and $\varepsilon=0$.

Now suppose $H$ is coupled to a continuous variable system such that the resulting interaction Hamiltonian is $H'=H\otimes \hat{z}$. Suppose the first register is prepared in some initial state $\ket{\psi_0}$ and the ancilla system is prepared in the continuous-variable quantum state

\ket{\bar{0}}_c=\int_{a}^{b}dz~\sqrt{\dfrac{c(z)}{\norm{c}_1}}\ket{z},

where $\norm{c}_1=\int_{a}^b dz~ |c(z)|$. For instance, $\hat{z}$ can represent a degree of freedom (position or momentum) of a one-dimensional quantum Harmonic oscillator, and the state $\ket{\bar{0}}_c$, could be its ground state (a Gaussian), a free resource state for continuous variable systems. For several of our applications, we shall see that this is indeed the case.

Now we shall simply evolve the system according to the interaction Hamiltonian $H'$ to obtain

\begin{align}
\ket{\eta_t}=e^{-i\tilde{H}t}\ket{\psi_0}\ket{\bar{0}}_c&=\int_{a}^{b} dz~\sqrt{\dfrac{c(z)}{\norm{c}_1}} e^{-iHtz}\ket{\psi_0}\ket{z}\\
&=\dfrac{1}{\norm{c}_1}\int_{a}^{b} dz~c(z) e^{-iHtz}\ket{\psi_0}\ket{\bar{0}}_c+\ket{\Phi}^\perp,
\label{eq:final-state-analog-lcu-1}
\end{align}

where $\ket{\Phi}^{\perp}$ is a quantum state (not normalized) such that $(I\otimes \ket{\bar{0}}_c\bra{\bar{0}}_c)\ket{\Phi}^{\perp}=0$. Thus, we have prepared a quantum state that is $O(\varepsilon/\norm{c}_1)$-close to

\begin{equation}
\ket{\psi}=\dfrac{f(H)}{\norm{c}_1}\ket{\psi_0}\ket{\bar{0}}_c+\ket{\Phi}^\perp.
\end{equation}

Now post-selecting on having $\ket{0}_c$ in the second register we obtain  a state that is $\varepsilon$-close to $f(H)\ket{\psi_0}/\norm{f(H)\ket{\psi_0}}$ in the first register with probability $\norm{f(H)\ket{\psi_0}}^2/\norm{c}^2_1$. We will use this procedure to develop an analog quantum algorithm for preparing ground states of Hamiltonians in Sec. <ref>.

This continuous-time algorithm can be naturally generalized to the scenario where we want to implement $f(H)$ for some function $f:[-1,1]\mapsto \mathbb{R}$ such that

\left|f(x)-\int_{a_1}^{b_1}dz_1~c(z_1)\int_{a_2}^{b_2}dz_2~c(z_2)\cdots\int_{a_k}^{b_k}dz_k~c(z_k)e^{-itxz_1z_2\cdots z_k}\right|\leq \varepsilon.

This can be implemented by coupling the Hamiltonian $H$ with $k$ different ancillary continuous-variable systems such that the effective interaction Hamiltonian is $\tilde{H}=H\otimes \hat{z}_1\otimes\cdots\otimes\hat{z}_k$. The $j$-th ancilla system is prepared in the quantum state

\ket{\bar{0}}_{c_j}=\int_{a_j}^{b_j} dz_j~\sqrt{\dfrac{c(z_j)}{\norm{c_j}_1}}\ket{z_j}.

Then by evolving the initial state $\ket{\psi_0}\ket{\bar{0}}_{c_1}\cdots\ket{\bar{0}}_{c_k}$ according to $\tilde{H}$ for time $t$ results in a quantum state that is $O\left(\frac{\varepsilon}{\Pi_{j=1}^k\norm{c_j}_1}\right)$-close to

\begin{equation}
\label{eq:final-state-analog-lcu-2}
\ket{\eta_t}=\dfrac{f(H)}{\Pi_{j=1}^k\norm{c_j}_1}\ket{\psi_0}\ket{\bar{0}}_{c_1}\cdots\ket{\bar{0}}_{c_k}+\ket{\Phi}^\perp.
\end{equation}

In Sec. <ref>, our analog quantum linear systems algorithm requires coupling the system Hamiltonian to two ancillary systems, which is captured by this generalization of analog LCU. Interestingly, for the applications we consider, the ancillary states are the ground or the first excited state of a one-dimensional quantum Harmonic oscillator or the ground state of a “particle in a ring”.

§.§ Robustness of expectation values of observables

In this section, we develop general results on the robustness of expectation values of observables which we shall use for both the “Single-Ancilla LCU" (Sec. <ref>) and the “Ancilla-free LCU” (Sec. <ref>)  approaches.

Consider that there exist two operators $P$ and $Q$ such that $\norm{P-Q}\leq \varepsilon$. In this section, we demonstrate that the expectation value of $O$ with respect to $P\rho P^{\dag}$ is not far off from the expectation value of $O$ with respect to $Q\rho Q^{\dag}$, for any density matrix $\rho$. More precisely, we prove

\left|\Tr[OP\rho P^{\dag}]-\Tr[OQ\rho Q^{\dag}]\right|\leq 3\norm{P}\norm{O}\varepsilon.

In order to prove this result, we need to use the tracial version Hölder's inequality which is stated below for completeness:

Define two operators $A$ and $B$ and parameters $p,q\in [1,\infty)$ such that $1/p+1/q =1 $. Then the following holds:

\Tr[A^{\dag}B]\leq \norm{A}_p \norm{B}_q.

Now we are in a position to formally state the main result of this section.

Suppose $P$ and $Q$ are operators such that $\norm{P-Q}\leq \varepsilon$ for some $\varepsilon\in [0,1]$. Furthermore, let $\rho$ be any density matrix and $O$ be some Hermitian operator with spectral norm $\norm{O}$. Then, if $\norm{P}\geq 1$, the following holds:

\left|\Tr[OP\rho P^{\dag}]-\Tr[OQ\rho Q^{\dag}]\right| \leq 3\norm{O}\norm{P}\varepsilon.

For the operators $P$ and $Q$, we have

\begin{align}
(Q-P)\rho(P^{\dag}-Q^{\dag})=Q\rho (P^{\dag}-Q^{\dag})-P\rho(P^{\dag}-Q^{\dag})
\end{align}

Now adding and subtracting $P\rho P^{\dag}$ in the RHS we obtain

\begin{align}
(Q-P)\rho(P^{\dag}-Q^{\dag})+P\rho(P^{\dag}-Q^{\dag})&= Q\rho (P^{\dag}-Q^{\dag})+P\rho P^{\dag}-P\rho P^{\dag}\\
&= P\rho P^{\dag} - Q\rho Q^{\dag} - (P-Q)\rho P^{\dag}
\end{align}

This gives us

\begin{equation}
P\rho P^{\dag} - Q\rho Q^{\dag} = (Q-P)\rho(P^{\dag}-Q^{\dag})+P\rho(P^{\dag}-Q^{\dag})+(P-Q)\rho P^{\dag}.
\end{equation}

We multiply $O$ to the left of each term on both sides of this equation and then take trace. Thus, we have

\begin{align}
\label{eq:distance-operators-restated}
\left|\Tr[O P\rho P^{\dag}] - \Tr[O Q\rho Q^{\dag}]\right| &= \left|\Tr[O (Q-P)\rho(P^{\dag}-Q^{\dag})]+\Tr[O P\rho(P^{\dag}- Q^{\dag})]+\Tr[(P-Q)\rho P^{\dag}]\right|\\
&\leq \left|\Tr[O (Q-P)\rho(P^{\dag}-Q^{\dag})]\right|+\left|\Tr[O P\rho(P^{\dag}- Q^{\dag})]\right|+\left|\Tr[(P-Q)\rho P^{\dag}]\right|,
\end{align}

Now for each term in the RHS, we invoke Theorem <ref> with $p=q=2$ to obtain

\begin{align}
\left|\Tr[O P\rho P^{\dag}] - \Tr[O Q\rho Q^{\dag}]\right| &\leq \norm{O}\norm{P-Q}^2+2\norm{O}\norm{P}\norm{P-Q}~~~~~~~~~[~\mathrm{As~}\norm{\rho}=1]\\
&\leq \varepsilon^2\norm{O}+2\norm{O}\norm{P}\varepsilon~~~~~~~~~~~~~~~~~~~~~~~~~~~~[~\mathrm{As~} \norm{P-Q}\leq \varepsilon]\\
&\leq 3\varepsilon\norm{O}\norm{P}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~[~\mathrm{As~} \norm{P}\geq 1]
\end{align}

It is easy to see why Theorem <ref> is useful to develop robust versions of the new approaches to LCU.  Typically, $f(H)$ is often not exactly equal to a linear combination of unitaries but is $\varepsilon$-close to it. Formally, $\norm{f(H)-f(H')}\leq \varepsilon$, where $f(H')$ can be exactly expressed as an LCU. Consequently, for the variants of LCU that we develop in the subsequent sections, we will be estimating $\Tr[O f(H')\rho f(H')^{\dag}]$. But, by Theorem <ref>, we can always bound

\left|\Tr[O f(H)\rho f(H)^{\dag}]-\Tr[O f(H')\rho f(H')^{\dag}]\right|\leq 3\norm{O}\norm{f(H)}\varepsilon.

We shall be using this result in the subsequent sections.

§.§ Single-Ancilla LCU: Controlled application of randomly sampled unitaries followed by single-shot measurement

In this section, we describe the “Single-Ancilla LCU" technique, which allows us to sample from quantum states obtained by applying LCU. This approach is a generalization of the Hamiltonian simulation algorithm of [9].

Suppose we are given a Hermitian matrix $H$ and we wish to implement $f(H)$, which can be approximated by a linear combination of unitaries. That is, for some $\gamma\in [0,1)$,

\norm{f(H)-\sum_{j=1}^{M} c_j U_j}\leq \gamma,

for unitaries $U_j$ and $c_j\in\mathbb{R}$. Let us define the $\ell_1$-norm of the LCU coefficients as $\norm{c}_1=\sum_{j=1}^{M} c_j$. Furthermore, suppose we have access to the quantum circuit for $U_j$. Then given any initial state $\ket{\psi_0}$, by the addition of just one ancilla qubit, we can estimate the expectation values $\Tr[O f(H)\ket{\psi_0}\bra{\psi_0} f(H)^{\dag}]$ up to arbitrary accuracy, for any observable $O$. This is a randomized quantum algorithm that is applicable to early fault-tolerant quantum devices: we use a short-depth quantum circuit, from which we sample repeatedly. We depict the resulting quantum circuit in Fig. <ref> and state the algorithm formally in Algorithm <ref>. The overall procedure is simple: The ancilla qubit is prepared in the state $\ket{+}$ so that the overall initial state is $\rho=\ket{+}\bra{+}\otimes \ket{\psi_0}\bra{\psi_0}$. We sample two unitaries $V_1$ and $V_2$ according to $\{U_j, c_j/\norm{c}_1\}$ and then implement controlled and anti-controlled versions of $V_1$ and $V_2$, respectively. Finally, we make a measurement of the observable $X\otimes O$ and store the outcome.

$\texttt{Single-Ancilla-LCU~}(O, \{c_j,U_j\},\ket{\psi_0}, T$)

1.  Prepare the state $\rho_0=\ket{+}\bra{+}\otimes \ket{\psi_0}\bra{\psi_0}$.

2.  Obtain i.i.d. samples $V_1, V_2$ from the distribution $\left\{U_j, \dfrac{c_j}{\norm{c}}\right\}$.

3.  For $\tilde{V}_1=\ket{0}\bra{0}\otimes I+\ket{1}\bra{1}\otimes V_1$ and $\tilde{V}_2=\ket{0}\bra{0}\otimes V_2\ket{1}\bra{1}\otimes I$, perform measurement of
$(X\otimes O)$ on the state

\rho'= \tilde{V}_1 \tilde{V}_2\rho_0 \tilde{V}^{\dag}_2\tilde{V}^{\dag}_1.

4.  For the $j^{\mathrm{th}}$ iteration, store into $\mu_j$, the outcome of the measurement in Step 4.

5.  Repeat Steps 1, 2, 3, and 4, a total of $T$ times.

6.  Output

\mu=\dfrac{\norm{c}^2_1}{T}\sum_{j=1}^{T} \mu_j.

We then sample from this quantum circuit $T$ times and calculate the mean of all the outcomes. This is the final output of Algorithm <ref>. We show formally via Theorem <ref>, that if the number of samples obtained $T$ is large enough, the sample mean of the outcomes converges to the expectation value we intend to estimate. Formally we define the following theorem,

Let $f(H')=\sum_{j}c_j U_j$. First observe from Algorithm <ref> that the initial state $\rho$ transforms to

\begin{align}
\rho'&= \tilde{V}_1 \tilde{V}_2 \rho \tilde{V}^{\dag}_2\tilde{V}^{\dag}_1\\
&=\dfrac{1}{2}\left[\ket{0}\bra{0}\otimes V_2 \rho_0 V^{\dag}_2 +\ket{0}\bra{1}\otimes V_2 \rho_0 V^{\dag}_1+\ket{1}\bra{0}\otimes V_1 \rho_0 V^{\dag}_2+\ket{1}\bra{1}\otimes V_1 \rho_0 V^{\dag}_1\right].
\end{align}

So after measuring the observable $X\otimes O$, we have

\Tr\left[(X\otimes O)\rho'\right]=\dfrac{1}{2}\left[O\left(V_1 \rho_0 V^{\dag}_2+ V_2 \rho_0 V^{\dag}_1\right)\right].

The expected values,

\mathbb{E}\left[V_1\right]=\mathbb{E}\left[V_2\right]=\dfrac{1}{\norm{c}_1}\sum_{j}c_j U_j.

So, the expected outcome of the $j^{\mathrm{th}}$ iteration is

\mathbb{E}\left[\mu_j\right]=\mathbb{E}\left[\Tr[(X\otimes O)\rho']\right]=\dfrac{1}{\norm{c}_1^2}\Tr[O~f(H')\rho_0 f(H')^{\dag}].

Next, we need to estimate two things:

(a) How fast does the sample mean $\mu=\sum_{j}\norm{c}^2_1\mu_j/T$ converge to its expectation value? For this, we use Hoeffding's inequality.

(b) What is the accuracy of the observation with respect to $f(H)$ as a function of the distance between $f(H)$ and $f(H')$? For this, we invoke Theorem <ref>.

Furthermore, observe that the POVM measurement yields some outcome of $O$ in the range $[-\norm{O},\norm{O}]$. So each random variable lies in the range

-{\norm{O}\norm{c}^2_1 \leq \norm{c}^2_1\mu_j\leq -\norm{O}\norm{c}^2_1}.

We evaluate (a), by using Hoeffding's inequality. We obtain

\begin{align}
Pr\left[\left|\mu - \Tr[O~f(H')\rho_0 f(H')^{\dag}]\right|\geq \varepsilon/2\right] \leq 2\exp\left[-\dfrac{T\varepsilon^2}{8\norm{c}^4_1\norm{O}^2}\right].
\end{align}

This immediately gives us that for

\begin{equation}
\label{eq:number-of-repititions}
T\geq \dfrac{8\norm{O}^2\ln(2/\delta)\norm{c}_1^4}{\varepsilon^2},
\end{equation}

\left|\mu -\ Tr[O~f(H')\rho_0 f(H')^{\dag}]\right| \leq \varepsilon/2,

with probability at least $1-\delta$. Now, in order to evaluate (b), we first apply triangle inequality, we obtain

\begin{align}
\left|\mu - \Tr[O~f(H)\rho_0 f(H)^{\dag}]\right| \leq &\left|\mu - \Tr[O~f(H')\rho_0 f(H')^{\dag}]\right|+ \\
&\left|\Tr[O~f(H)\rho_0 f(H)^{\dag}]-
\Tr[O~f(H')\rho_0 f(H')^{\dag}]\right|.
\end{align}

The first term in the RHS of the above inequality is upper bounded by $\varepsilon/2$. In order to bound the second term, note that $\norm{f(H)-f(H')}\leq \gamma$. For any such operators that are at most $\gamma$-separated, we can use Theorem <ref> to obtain:

\left|\Tr[O~f(H)\rho_0 f(H)^{\dag}]-
\Tr[O~f(H')\rho_0 f(H')^{\dag}]\right|\leq 3\norm{O}\norm{f(H)}\gamma\leq \varepsilon/2.

So, overall we have

\left|\mu - \Tr[O~f(H)\rho_0 f(H)^{\dag}]\right| \leq \varepsilon,

which completes the proof.

What is this procedure useful for?

In Sec. <ref> and Sec. <ref>, we show that this approach is useful to estimate expectation values $\Tr[~O f(H)\rho f(H)^{\dag}]$ for several functions of interest. In particular, we apply this to sample from the solution of quantum linear systems and also from the ground states of Hamiltonians. These quantum algorithms consume minimum quantum resources and rely on repeated classical sampling from a short-depth quantum circuit and hence, are applicable for early fault-tolerant quantum computers.

Next, we discuss the “Ancilla-free LCU” procedure.

§.§ Ancilla-free LCU: Randomized unitary sampling

As in the previous section, suppose we are given a Hermitian matrix $H$ and we wish to implement $f(H)$, which can be approximated by a linear combination of unitaries. That is, for some $\gamma\in [0,1)$,

\norm{f(H)-\sum_{j=1}^{M} c_j U_j}\leq \gamma,

for unitaries $U_j$, $c_j\in\mathbb{R}$, and $\norm{c}_1=\sum_{j=1}^{M} c_j$. In this section, we formally prove that if we are interested in the projection of $f(H)\ket{\psi_0}$ in some subspace of interest, then it suffices to sample $U_j$ according to the distribution of the LCU coefficients, i.e. $\mathcal{D}\sim \{c_j/\norm{c}_1\}$. This is the key idea behind the “Ancilla-free LCU” technique.

For some initial state $\rho_0=\ket{\psi_0}\bra{\psi_0}$, this procedure prepares the average density matrix

\rho=\dfrac{1}{\norm{c}_1}\sum_{j=1}^{M} c_j U_j \rho_0 U^{\dag}_j.

Then, for some projector $\Pi$ onto the subspace of interest, the projection of $\rho$ is at least as large as the projection of $f(H)\ket{\psi_0}$. We  formally prove this via the following theorem:

Let $f(H')=\sum_{j=1}^{M} c_j U_j$. Then, the standard LCU procedure would have implemented the quantum state

\ket{\psi_t}=\ket{\bar{0}} \dfrac{f(H')}{\norm{c}_1}\ket{\psi_0} +\ket{\Phi^{\perp}}.

Then, we have

\begin{align}
\Tr[\Pi\rho]&=\Tr[(I\otimes \Pi)\ket{\psi_t}\bra{\psi_t}]\\
&=\bra{\psi_t}\left(\ket{\bar{0}}\bra{\bar{0}}\otimes \Pi\right)\ket{\psi_t}+\bra{\psi_t}\left[(I-\ket{\bar{0}}\bra{\bar{0}})\otimes \Pi\right]\ket{\psi_t}\\
&= \dfrac{1}{\norm{c}^2_1}\bra{\psi_0}f(H')^{\dag} \Pi f(H')\ket{\psi_0}+\bra{\Phi^\perp}\left[(I-\ket{\bar{0}}\bra{\bar{0}})\otimes\Pi \right]\ket{\Phi^\perp}\\
&\geq \dfrac{1}{\norm{c}^2_1}\bra{\psi_0}f(H')^{\dag} \Pi f(H')\ket{\psi_0}\\
&\geq \dfrac{1}{\norm{c}^2_1}\left(\bra{\psi_0}f(H)^{\dag} \Pi f(H)\ket{\psi_0}-\varepsilon\right),
\end{align}

where in the last line we have invoked Theorem <ref>. This completes the proof.

A similar result naturally extends to the analog LCU framework as well. Thus, replacing the LCU procedure by random unitary sampling suffices in cases where we are interested in the projection of $f(H)\ket{\psi_0}$ in some subspace, i.e. say we wish to
measure this state in some basis. Then we may instead measure $\rho$, for which the probabilities would be at least as large.

This is typically the case for several quantum walk-based problems: the goal is to prepare a quantum state that has a good overlap with the marked nodes (some subset of the nodes) of the underlying graph. This problem can be tackled by discrete-time quantum walks as well as by continuous-time quantum walks.

Briefly, $f(D)=D^t$ represents $t$-steps of a discrete-time random walk on $D$, the discriminant matrix of some (reversible) Markov chain $P$. Let $\Pi_M$ be a projector onto the space spanned by the marked nodes of $P$. Then, one can express $f(D)=D^{T}$ as a linear combination of discrete-time quantum walk steps for some $T=O(HT)$, which has roughly $\sqrt{T}$ terms ($HT$ is the hitting time of the discrete-time random walk on $D$). Then, instead of implementing the complete LCU procedure, one can apply $k$-steps of a discrete-time quantum walk, where $k$ is sampled according to the distribution of the LCU coefficients. This results in $\rho$ such that $\Tr[\Pi_M\rho]$, the probability of observing a marked vertex after $\sqrt{T}$ quantum walk steps is at least as large as the classical quantity $\braket{\psi_0|D^T \Pi_M D^T|\psi_0}$, which corresponds to the distribution after $T$-steps of a classical random walk. The classical quantity can be proven to be $Ω̃(1)$ for some specific $|ψ_0⟩$, using the framework of interpolated Markov chains (following the work of Amabainis et al [31]), which proves optimality of the quantum spatial search algorithm. Similar results can also be shown with respect to continuous-time random walks, wherein the LCU decomposition $f(D)=e^T(D-I)$ becomes relevant. For continuous-time quantum walks [32], on the other hand, $f(D)=e^T(D^2-I)$. Thus, the ``Ancilla-Free'' LCU approach is the broad framework to tackle quantum walk-based algorithms. We shall discuss these ideas in Sec.~\ref{sec:quantum-random-walks}
\section{Applications to Ground state preparation}
\label{sec:gsp}
The ground state preparation (GSP) problem can be stated as follows: given access to a Hamiltonian and an initial state, find a quantum algorithm that outputs the ground state of the Hamiltonian with high precision. Generally, this problem is known to be computationally hard, even for a quantum computer [68]. However, it finds widespread applications across physics and computer science.  As a result, novel quantum algorithms for GSP that improve upon existing ones are of extreme importance and interest. In this section, we will develop quantum algorithms for preparing ground states of Hamiltonians as well for sampling from the ground state. We will find that both the ``Analog LCU'' and the ``Single-Ancilla LCU'' approaches, introduced in Sec.~\ref{sec:new-approaches-LCU}, are useful for this problem.

This section is organized as follows. We begin by formally describing the GSP problem. Next, we develop an analog quantum algorithm for GSP using the ``Analog LCU'' framework. Then, we use the ``Single-Ancilla LCU'' technique to sample from the ground state. Finally, as a bonus, we provide a quantum algorithm for GSP for fully fault-tolerant quantum computers using the framework of QSVT.

We start by describing the ground state preparation problem.
\\~\\
\textbf{The ground state Preparation problem:~} The set up of the problem is similar to prior works [34, 35, 32]. Suppose we have a Hamiltonian $H$ with ground state $|v_0⟩$ and ground energy $λ_0$, and assume that we are given a lower bound on the gap between the ground state and the first excited state of $H$ (spectral gap), i.e.\ we have knowledge of $Δ$ such that $|λ_1-λ_0|≥Δ$.\\
For clarity of exposition, we assume that the ground space of $H$ is non-degenerate. If this is not the case, e.g.~if the degeneracy of the ground space is $d$ and is spanned by mutually orthonormal eigenstates ${|v^(ℓ)_0⟩}_l=1^d$, then we will be preparing a quantum state $|v_0⟩$ which is a projection onto the ground space given by
\ket{v_0}=\dfrac{1}{\sqrt{\sum_{\ell=1}^{d} |c^{(\ell)}_0|^2}}\sum_{\ell=1}^d c^{(\ell)}_0\ket{v^{(\ell)}_0}.
In addition, suppose we have access to some initial state $|ψ_0⟩$ and a lower bound on the overlap $|⟨ψ_0|v_0||⟩=c_0≥η$. \\
Furthermore, for some desired accuracy $∈(0,1)$, we will assume that we know the value of the ground energy to some precision parameter $_g$ such that $_g = (Δ/√(log1/η))$. That is, we know some $E_0$ such that
\begin{equation}
% \label{eqsup:ground-energy-range}
|\lambda_0-E_0| \leq \eps_g.
\end{equation}
By implementing $H-(E_0-_g)I$, we ensure that $0≤λ_0≤2_g$. This transformation also ensures that the lower bound for the spectral gap of $H$ remains $Δ$. If also an upper bound on the maximum eigenvalue of $H$ is known, then we can actually assume that the spectrum of $H$ is in $[0,1]$.
\\~\\
\subsection{Applying Analog LCU: A continuous-time quantum algorithm for ground state preparation}
\label{subsec:analog-gsp}
In this section, we will use the ``Analog LCU'' framework to develop an analog quantum algorithm for the GSP problem. This algorithm was described in the Supplemental Material of [32]. Here, we place it in the broader context of the ``Analog LCU'' framework. Moreover, it will serve as useful intuition for the ``Ancilla free LCU" and the ``Single-Ancilla LCU" approaches to the GSP problem, which we discuss in the subsequent sections.

Consider some quantum system in state $|ψ_0⟩$ coupled to an ancillary system in a Gaussian state
\begin{equation}
\label{eq:gaussian-state}
\ket{\psi_g}=\int_{-\infty}^{+\infty} \dfrac{dz}{(2\pi)^{1/4}} e^{-z^2/4}\ket{z}.
\end{equation}
The Gaussian state is typically easy to prepare in this setting. This state can be seen as the ground state of a one-dimensional quantum harmonic oscillator. The coupling is done via interaction Hamiltonian $H'=H⊗ẑ$, where $ẑ$ corresponds to the position (or momentum) operator. Evolving $|ψ_0⟩|ψ_g⟩$, under $H'$ for a time $t$ results in the state
\begin{align}
\ket{\eta_t}&=e^{-itH'}\ket{\psi_0}\ket{\psi_g} \nonumber \\
\label{eqmain:evolved-state-fast-forwarding-1}
&=\int_{-\infty}^{+\infty}\dfrac{dz}{(2\pi)^{1/4}}e^{-z^2/4}e^{-itHz}\ket{\psi_0}\ket{z} \\
&=\int_{-\infty}^{+\infty}\dfrac{dz}{\sqrt{2\pi}}e^{-z^2/2}e^{-itHz}\ket{\psi_0}\ket{\psi_g}+\ket{\Phi}^\perp, \nonumber
\end{align}
where $|Φ⟩^⊥$ is a quantum state with the ancillary system being orthogonal to $|ψ_g⟩$. Now the Fourier transform of a Gaussian is a Gaussian, i.e.\ we have for any $x ∈ℝ$,
\begin{equation}
\label{eq:hubbard-stratonovich}
e^{-y^2/2}=\int_{-\infty}^{\infty} \dfrac{dz}{\sqrt{2\pi}}~e^{-z^2/2}e^{-iyz}.
\end{equation}
So using this we obtain
\begin{equation}
\label{eqmain:evolved-state-fast-forwarding-3}
\ket{\eta_t}=e^{-t^2H^2/2}\ket{\psi_0}\ket{\psi_g}+\ket{\Phi}^\perp.
\end{equation}
By post-selecting on obtaining $|ψ_g⟩$ in the second register, we are able to prepare a quantum state proportional to $e^-t^2H^2/2|ψ_0⟩$ in the first register. Now we formally state the ground state preparation algorithm and analyze its complexity, via the following lemma.
\begin{lemma}
\label{lem:analog-gsp-ham}
Suppose $\eps\in (0,1)$ and $\eta\in (0, 1/\sqrt{2}]$. Furthermore, suppose we have a Hamiltonian $H$ with ground state $\ket{v_0}$ with $\Delta$ being a lower bound on the spectral gap. Also, the ground state energy of $H$ is known up to a precision $\eps_g\in \Oo\left(\Delta/\sqrt{\log\frac{1}{\eta\eps}}\right)$. Then, given an initial state $\ket{\psi_0}$ satisfying $|\braket{\psi_0|v_0}|\geq \eta$, we output, with probability $\Oo(\eta^2)$, a state $\ket{\phi}$ such that $\nrm{\ket{\phi}-\ket{v_0}}\leq \eps$ by evolving the Hamiltonian $H'=H\otimes \hat{z}$ for time
T = \Oo\left(\dfrac{1}{\Delta}\sqrt{\log \left(\dfrac{1}{\eta\eps}\right)}\right).
\end{lemma}
\begin{proof}
We shift the overall eigenvalues of $H$ by $-E_0+\eps_g$ as explained previously. So now suppose $H$ has a spectral decomposition is $H=\sum_{j}\lambda_j\ket{v_j}\bra{v_j}$ with eigenvalues $\lambda_j\in (0,1]$ and particularly $0\leq \lambda_0\leq 2\eps_g$.

Then the input quantum state, when expressed in the eigenbasis of $H$ is written as $\ket{\psi_0}=\sum_{j}c_j\ket{v_j}$. Without loss of generality, we assume that $c_0$ is real and positive and follow the analog procedure outlined in the main article. In fact, from Eq.~(3) of the letter, we obtain that
\begin{equation}
\ket{\eta_t}=e^{-tH^2}\ket{\psi_0}\ket{\psi_g}+\ket{\Phi}^\perp. \nonumber
\end{equation}
By post-selecting on obtaining the Gaussian state $\ket{\psi_g}$ on the second register, we are left with the following quantum state in the first register after a time $\sqrt{2t}$
\begin{equation}
% \label{eqsup:ground-state-prepare}
\ket{\phi}=\dfrac{e^{-tH^2}\ket{\psi_0}}{\sqrt{\braket{\psi_0|e^{-2tH^2}|\psi_0}}},
\end{equation}
with probability $\| e^{-tH^2}\ket{\psi_0}\ket{\psi_g} \|^2 = \braket{\psi_0|e^{-2tH^2}|\psi_0}$. Expressing $\ket{\phi}$ in the eigenbasis of $H$, we obtain
\begin{equation}
% \label{eqsup:phi-in-ebasis-of-H}
\ket{\phi}=\dfrac{c_0e^{-t\lambda^2_0}}{\sqrt{\braket{\psi_0|e^{-2tH^2}|\psi_0}}}\left[\ket{v_0}+\sum_{j\geq 1}\dfrac{c_j}{c_0} e^{-t(\lambda^2_j-\lambda^2_0)}\ket{v_j}\right],
\end{equation}
where the normalization factor
\begin{equation}
% \label{eqsup:succ-prob-gsp}
\braket{\psi_0|e^{-2tH^2}|\psi_0}=|c_0|^2 e^{-2t\lambda^2_0}\left[1+\sum_{j\geq 1}\frac{|c_j|^2}{|c_0|^2}e^{-2t(\lambda^2_j-\lambda^2_0)}\right].
\end{equation}
Now, we intend to choose a value of $t$, so that $\ket{\phi}$ and the ground state $\ket{v_0}$ are $\eps$-close to each other in $\ell_2$-norm. We have,
\begin{equation}
\nrm{\ket{\phi}-\ket{v_0}}^2 = 2-2\braket{\phi|v_0},
\end{equation}
where we use the fact that $c_0\geq 0$ which implies $\braket{\phi|v_0}> 0$. Thus we have,
\begin{align}
\braket{\phi|v_0}&=\left[1+\sum_{j\geq 1} \dfrac{|c_j|^2}{|c_0|^2}e^{-2t(\lambda_j^2-\lambda_0^2)}\right]^{-1/2} \nonumber\\
&\geq \left[1+\dfrac{1-\eta^2}{\eta^2} e^{-2t(\lambda_1^2-\lambda_0^2)}\right]^{-1/2}\\
&\geq 1-\dfrac{1-\eta^2}{2\eta^2} e^{-2t(\lambda_1^2-\lambda_0^2)},
\label{eq:gstate-phi-overlap}
\end{align}
which gives
\begin{equation}
\label{eq:distance-ground-state}
\nrm{\ket{\phi}-\ket{v_0}}^2\leq \dfrac{1-\eta^2}{\eta^2} e^{-2t(\lambda_1^2-\lambda_0^2)}\leq \dfrac{1-\eta^2}{\eta^2} e^{-2t(\lambda_1-\lambda_0)^2}.
\end{equation}
Thus by choosing any value of $t$ such that
\begin{align}
\end{align}
we ensure that $\nrm{\ket{\phi}-\ket{v_0}}\leq \eps$. The total evolution time of the interaction Hamiltonian is $T=\sqrt{2t}$, which is
\begin{equation}
% \label{eqsup:running-time-gsp}
T= \Oo\left(\dfrac{1}{\Delta}\sqrt{\log \left(\dfrac{1}{\eta\eps}\right)}\right).
\end{equation}
Now we have,
\begin{equation}
\braket{\psi_0|e^{-2tH^2}|\psi_0}= \dfrac{|c_0|^2 e^{-2t\lambda^2_0}}{|\braket{\phi|\psi_0}|^2}\geq |c_0|^2 e^{-2t\lambda^2_0}.
\end{equation}
So the success probability of our algorithm
\begin{equation}
\label{eq:succ-prob-gsp-2}
\braket{\psi_0|e^{-T^2H^2}|\psi_0}\geq |c_0|^2 e^{-T^2\lambda^2_0} \geq \Oo(\eta^2),
\end{equation}
where we have used the fact that $\lambda_0\leq 2\eps_g$ and so $T\lambda_0 = \Oo(1)$.
\end{proof}
Overall, this physically motivated quantum algorithm is significantly simpler than implementing standard LCU in the circuit model. Moreover, hybrid qubit-qumode systems are currently being engineered in a number of quantum technological platforms. In the future, we intend to provide an experimental proposal to implement ``Analog LCU" on experimental platforms such as ion traps or superconducting systems.

Next, we move on to applications of the other approaches of LCU for the GSP problem. Next, we describe how the ``Single-Ancilla LCU'' technique can be used to develop a randomized quantum algorithm for sampling from ground states.
\subsection{Applying Single-Ancilla LCU: Sampling from the ground states of Hamiltonians}
\label{subsec:gsp-single-ancilla}
In this section, suppose that $H\in \mathbb{R}^{N\times N}$ such that the ground state of $H$ is $\ket{v_0}$. We intend to use the techniques of Sec.~\ref{subsec:one-ancilla-LCU}, in particular Algorithm \ref{algo:randomized-time-evolution} and Theorem \ref{thm:randomized-time-evolution} to estimate the quantity $\braket{v_0|O|v_0}$, upto $\varepsilon$-accuracy for some observable $O$. As promised, our quantum algorithm will require sampling from a short-depth quantum circuit of Fig.~\ref{fig:single-ancilla-circuit} and hence, will be applicable to early fault-tolerant quantum computers.

The basic idea is to take advantage of the fact the LCU decomposition of $f(H)=e^{-tH^2}$. Following the analog quantum algorithm in Sec.~\ref{subsec:analog-LCU}, we already know that $f(H)\ket{\psi_0}$ helps in preparing the ground state of $H$, for an appropriate choice of $t$. In this section, we consider a discretized version of this LCU decomposition, i.e.\ we express $f(H)\approx \sum_{j=1}^{T} c_j e^{-ijH}$, where $T$ scales as roughly $\sqrt{t}$.  This allows us to use Algorithm \ref{algo:randomized-time-evolution} to estimate $\braket{v_0|O|v_0}$.

Let us begin by obtaining an $\gamma$-accurate LCU decomposition of $f(H)=e^{-tH^2}$. This decomposition has already shown up in prior works [36, 66]. We formally state this via the following Lemma which we prove, for completeness in the Appendix.

\begin{restatable}[LCU decomposition of $e^{-tH^2}$]{lemma}{lemmalcuexp}
\label{lem:lcu-decomp-exp}
Let $0<\gamma<1$ and consider a Hamiltonian $H$ of unit spectral norm. Furthermore, for any $t>1$, let us define
\begin{equation*}
X_M= \sum_{j=-M}^M c_j e^{-i j \delta_t \sqrt{2 t} H },
\end{equation*}
where $M=\left \lceil{\sqrt{2}\left(\sqrt{t}+\sqrt{\log(5/\gamma)}\right)\sqrt{\log(4/\gamma)}}\right \rceil $, $
\delta_t= \left(\sqrt{2t}+\sqrt{2\log(5/\gamma)}\right)^{-1}$ and,
c_j=\frac{\delta_t}{\sqrt{2 \pi}} e^{- j^2 \delta_t^2/2}.
\nrm{X_M-e^{-t H^2}} \leq \gamma.
\end{restatable}
\begin{proof}
We prove this in the Appendix (See Sec.~\ref{sec:appendix-proof-LCUdecomp})
\end{proof}

We shall use this lemma to develop our randomized quantum algorithm. First, observe that the $\ell_1$-norm of the LCU coefficients can be upper bounded by a constant. In fact,
\begin{align}
||c||_1&=\sum_{j=M}^M |c_j|\\&\leq |c_0|+2\sum_{j=1}^{\infty}\frac{\delta_t}{\sqrt{2\pi}}e^{-j^2\delta_t^2/2}\\
&\leq |c_0|+2 \int_0^\infty\frac{e^{-x^2/2}}{\sqrt{2\pi}}dx
=1+|c_0|\leq 1+\delta_t=O(1).
\label{eq:l1-norm-exp}
\end{align}
We will now use these results to estimate $\braket{v_0|O|v_0}$ by invoking Algorithm \ref{algo:randomized-time-evolution}. First of all, notice that each iteration of our randomized quantum algorithm requires sampling $V_1, V_2$ according to $\{U_j, c_j/\norm{c}_1\}$. From Lemma \ref{lem:lcu-decomp-exp}, we know that each $U_j=e^{-ij\delta_t\sqrt{2t}H}$. Thus, the cost of implementing each $U_j$ would be the cost of simulating $H$. For this, any Hamiltonian simulation technique can be used. Some of the recent randomized Trotter-based approaches [51, 38, 9, 40] will be more suitable for implementing this algorithm on early fault-tolerant quantum computers. For simplicity and clarity of exposition, we consider the Hamiltonian simulation subroutine as a black box and assume that  the cost of implementing the Hamiltonian simulation is upper bounded by the maximum time of evolution for $H$. From Lemma \ref{lem:lcu-decomp-exp}, the largest evolution time is given by
\begin{equation}
\label{eq:gstate-max-time-evolution}
\tau_{\max}=M \delta_t\sqrt{2t}=O\left(\sqrt{t\log\left(1/\gamma\right)}\right),
\end{equation}
where, for the GSP problem, the choice of $\sqrt{t}$ can be obtained from Lemma \ref{lem:analog-gsp-ham}.

So, for any observable $O$ and an initial state $\rho_0=\ket{\psi_0}\bra{\psi_0}$, each run of Algorithm \ref{algo:randomized-time-evolution} estimates some $\mu_j$, such that
\mathbb{E}\left[\mu_j\right]=\Tr[O~X_M\rho_0 X^{\dag}_M].
However, if $\rho_g=\ket{v_0}\bra{v_0}$ is the ground state of $H$, we want the output of Algorithm \ref{algo:randomized-time-evolution}, to be some $\mu$ such that
\left|\mu- Tr[O\rho_g]\right|\leq \varepsilon.
For this, the precision $\gamma$ must be judiciously chosen, using Theorem \ref{thm:distance-expectation} so that the output $\mu$ is $\varepsilon$-close to the desired expectation value. Furthermore, from Lemma \ref{lem:analog-gsp-ham}, we know that only the normalized quantum state
\dfrac{e^{-tH^2}\ket{\psi_0}}{\norm{e^{-tH^2}\ket{\psi_0}}},
is close to the ground state. Consequently, we need information about $\ell=\norm{e^{-tH^2}\ket{\psi_0}}$, which from Eq.~\eqref{eq:succ-prob-gsp-2}, we know to be at least $\eta$, the overlap between the ground state and the initial state.

For now, let us assume that we have exact knowledge of $\ell=\norm{e^{-tH^2}\ket{\psi_0}}$. We will relax this assumption later on. Then, we divide the output of $\mu$ of Algorithm \ref{algo:randomized-time-evolution} by $\ell^2$, i.e.\ obtain $\mu/\ell^2$. In fact, we prove via the following theorem that $\mu/\ell^2$ is a $\varepsilon$-accurate approximation of $\Tr[O\rho_g]$.
\begin{theorem}
\label{thm:ground-state-estimation}
Let $\varepsilon, \delta, \gamma \in (0,1)$ and $\eta \in (0,1/\sqrt{2}]$. Suppose $H$ is a Hermitian matrix with ground state $\ket{v_0}$ and let $\ket{\psi_0}$ be some initial state, prepared in cost $\tau_{\psi_0}$, such that $|\braket{v_0|\psi_0}|=\eta$. Let $O$ be some observable. Furthermore, for
define $\ell^2=\braket{\psi_0|e^{-2tH^2}|\psi_0}$ . Then if
\gamma \leq \dfrac{\varepsilon \eta^2}{12\norm{O}},
such that,
\norm{e^{-tH^2}-X_M}\leq \gamma,
T\geq O\left(\dfrac{\norm{O}^2\ln(2/\delta)}{\varepsilon^2\eta^4}\right)
then Algorithm \ref{algo:randomized-time-evolution} allows outputs, with probability at least $1-\delta$, a parameter $\mu$ such that
\left|\dfrac{\mu}{\ell^2}-\braket{v_0|O|v_0}\right|\leq \varepsilon,
using $T$ calls to the quantum circuit in Fig.~\ref{fig:single-ancilla-circuit}, where the cost of each such repetition is at most $\tau_{\max}+\tau_{\psi_0}$, where
\tau_{max}=O\left(\dfrac{1}{\Delta}\log\left(\dfrac{\norm{O}}{\varepsilon \eta}\right)\right).
\end{theorem}
\begin{proof}
We have
\begin{equation}
\norm{e^{-tH^2}-X_M}\leq \gamma,
\end{equation}
where the upper bound in $\gamma$ in the Theorem statement ensures $\gamma\leq \dfrac{\varepsilon \ell^2}{12\norm{O}}$ as $\ell>\eta$.

From Theorem \ref{thm:randomized-time-evolution}, we have that for this choice of $\gamma$,
\begin{equation}
\left|\dfrac{\mu}{\ell^2}-\dfrac{\Tr[O e^{-tH^2}\rho_0 e^{-tH^2}]}{\ell^2}\right|\leq \varepsilon/2.
\end{equation}
This implies that the sample complexity is
\begin{align}
\label{eq:sample-complexity-gstate}
T\geq O\left(\dfrac{\norm{O}^2\log(2/\delta)\norm{c}^4_1}{\varepsilon^2 l^4}\right) = O\left(\dfrac{\norm{O}^2\log(2/\delta)}{\varepsilon^2 \eta^4}\right),
\end{align}
where have used the fact that $\norm{c}_1=O(1)$ (from Eq.~\eqref{eq:l1-norm-exp}) and the bound $\ell\geq \eta$ (from Eq.~\eqref{eq:succ-prob-gsp-2}).

Now, Algorithm \ref{algo:randomized-time-evolution}, outputs $\mu$ such that
\begin{equation}
\dfrac{\mu}{\ell^2}=\dfrac{\norm{c}^2_1}{T \ell^2}\sum_{j=1}^{T} \mu_j.
\end{equation}
If $\rho_0=\ket{\psi_0}\bra{\psi_0}$ and $\rho_g=\ket{v_0}\bra{v_0}$, by triangle inequality,
\begin{align}
\left|\dfrac{\mu}{\ell^2}-\Tr[O\rho_g]\right|&\leq \left|\dfrac{\mu}{\ell^2}-\dfrac{\Tr[O e^{-tH^2}\rho_0 e^{-tH^2}]}{\ell^2}\right|+\left|\Tr[O\rho_g]-\dfrac{\Tr[O e^{-tH^2}\rho_0 e^{-tH^2}]}{\ell^2}\right|.
\label{eq:gstate-inequality-bound-distance}
\end{align}
We have already shown that the first expression in the RHS of Eq~\eqref{eq:gstate-inequality-bound-distance} is bounded by $\varepsilon/2$. In order to upper bound the second expression in the RHS of Eq.~\eqref{eq:gstate-inequality-bound-distance}, let us consider the normalized quantum state $\ket{\phi}$ as defined in Eq.~\eqref{eq:state-phi}, i.e.\
\ket{\phi}=\dfrac{e^{-tH^2}\ket{\psi_0}}{\norm{e^{-tH^2\ket{\psi_0}}}}
Then, we have:
\begin{align}
\left|\braket{v_0|O|v_0}-\braket{\phi|O|\phi}\right|&\leq \norm{O}_{\infty}\norm{\ket{\phi}\bra{\phi}-\ket{v_0}\bra{v_0}}_1~~~~[\text{~Using Lemma \ref{thm:holder} with~}p=\infty,~q=1]\\
&\leq 2\norm{O}\sqrt{1-|\braket{v_0|\phi}|}.
\label{eq:distance-expectation-gsp}
\end{align}
From Lemma \ref{lem:analog-gsp-ham}, we get that by choosing some
we ensure that
\begin{align}
\left|\braket{\phi|v_0}\right|\geq 1-\dfrac{\varepsilon^2}{16\norm{O}^2}.
\end{align}
Substituting this back in Eq.~\eqref{eq:distance-expectation-gsp} gives us,
\left|\braket{v_0|O|v_0}-\braket{\phi|O|\phi}\right|\leq \varepsilon/2,
as desired.
Finally, from the aforementioned choice of $t$, we also obtain $\tau_{\max}$, the maximum time of evolution required. This is now,
\begin{align}
\tau_{\max}=O\left(\sqrt{t\log\left(\dfrac{\norm{O}}{\varepsilon \eta^2}\right)}\right)=O\left(\dfrac{1}{\Delta}\log\left(\dfrac{\norm{O}}{\varepsilon \eta}\right)\right).
\end{align}
Combining the bounds we obtain for each of the expressions in the RHS of Eq.~\eqref{eq:gstate-inequality-bound-distance}, we finally obtain
\left|\dfrac{\mu}{\ell^2}-\Tr[O\rho_g]\right|\leq \varepsilon,
which completes the proof.
\end{proof}
In the previous theorem, we assumed exact knowledge of $\ell^2$. This requirement can be relaxed. In Sec.~\ref{sec:appendix-robustness-norm}, we prove that we only require the knowledge of an approximation $\tilde{\ell}$ of $\ell^2$ (Theorem \ref{thmsup:norm-approx}). We prove that if
\begin{equation}
\label{eq:approximation-norm}
\left|\tilde{\ell}-\ell^2\right|\leq \dfrac{\varepsilon\ell^2}{6\norm{O}},
\end{equation}
and the estimate $\mu$ in Theorem \ref{thm:ground-state-estimation} be such that
\left|\dfrac{\mu}{\ell^2} - \dfrac{\Tr[O f(H)\rho f(H)]}{\ell^2}\right|\leq \varepsilon/3,
then we have
\begin{equation}
\left|\dfrac{\mu}{\tilde{\ell}} - \dfrac{\Tr[O g(H)\rho g(H)]}{\ell^2}\right|\leq \varepsilon.
\end{equation}
Thus, we only need an approximation $\tilde{\ell}$ of $\ell^2$ that satisfies Eq.~\eqref{eq:approximation-norm}. Moreover, one can use the recently developed quantum algorithm in Ref.~[12] to obtain an approximation to $\ell^2$. In particular, the quantum algorithm in Ref.~[12] obtains samples from a Hadamard test circuit which is also applicable for early fault-tolerant quantum computers (short depth and needing only one ancilla qubit). So, one can make use of this algorithm to estimate $\tilde{\ell}$, an approximation to the normalization factor $\ell$, and then run our quantum algorithms. So, our assumption of the prior knowledge of $\tilde{\ell}$ is not a strong one. For our results, however, we assume that such an estimate $\tilde{\ell}$ has already been obtained. This gives us the main theorem of this section, which we state next:
\begin{theorem}[Sampling from ground states of Hamiltonians]
\label{thm:ground-state-sampling}
Let $\varepsilon, \delta, \gamma \in (0,1)$ and $\eta \in (0,1/\sqrt{2}]$. Suppose $H$ is a Hermitian matrix with ground state $\ket{v_0}$ and let $\ket{\psi_0}$ be some initial state, prepared in cost $\tau_{\psi_0}$, such that $|\braket{v_0|\psi_0}|=\eta$. Let $O$ be some observable. Also for,
define $\ell^2=\braket{\psi_0|e^{-2tH^2}|\psi_0}$ . Suppose we know $\tilde{\ell}$, an approximation of $\ell^2$ such that
\left|\tilde{\ell}-\ell^2\right|\leq \dfrac{\varepsilon\ell^2}{6\norm{O}}.
Then if
\gamma \leq \dfrac{\varepsilon \eta^2}{36\norm{O}},
such that,
\norm{e^{-tH^2}-X_M}\leq \gamma,
T\geq O\left(\dfrac{\norm{O}^2\ln(2/\delta)}{\varepsilon^2\eta^4}\right)
then Algorithm \ref{algo:randomized-time-evolution} allows outputs, with probability at least $1-\delta$, a parameter $\mu$ such that
\left|\dfrac{\mu}{\tilde{\ell}}-\braket{v_0|O|v_0}\right|\leq \varepsilon,
using $T$ calls to the quantum circuit in Fig.~\ref{fig:single-ancilla-circuit}, where the cost of each such repetition is at most $\tau_{\max}+\tau_{\psi_0}$, where
\tau_{max}=O\left(\dfrac{1}{\Delta}\log\left(\dfrac{\norm{O}}{\varepsilon \eta}\right)\right).
\end{theorem}
\begin{proof}
The upper bound on $\gamma$ ensures that
\gamma \leq \dfrac{\varepsilon \ell^2}{36\norm{O}},
as $\ell>\eta$. So from Theorem \ref{thm:randomized-time-evolution}, we have
\left|\dfrac{\mu}{\ell^2}-\dfrac{\Tr[O e^{-tH^2}\rho_0 e^{-tH^2}]}{\ell^2}\right|\leq \varepsilon/6.
Thus, combining this bound with Theorem \ref{thmsup:norm-approx}, it follows that as long as
\left|\tilde{\ell}-\ell^2\right|\leq \dfrac{\varepsilon\ell^2}{6\norm{O}},
we have
\left|\dfrac{\mu}{\tilde{\ell}}-\dfrac{\Tr[O e^{-tH^2}\rho_0 e^{-tH^2}]}{\ell^2}\right|\leq \varepsilon/2.
The rest of the proof follows exactly like the derivation of Theorem \ref{thm:ground-state-estimation}.
\end{proof}
\subsection{Ground state preparation using QSVT on fully fault-tolerant quantum computers}
\label{subsec:gsp-qsvt}
In this section, we provide a quantum algorithm for the GSP problem for fully fault-tolerant quantum computers. The key idea is to implement the function $f(H)=e^{-tH^2}$ in the circuit model. A straightforward approach would be to use the decomposition of $f(H)$ in Lemma \ref{lem:lcu-decomp-exp} and implement a standard LCU procedure. However, a more efficient approach would be to implement a polynomial approximation of $f(H)$ using QSVT. This is what we implement next.

A polynomial approximation of $e^{-tx^2}$ can be obtained for $x\in [-1,1]$ by modifying the polynomial in Lemma \ref{lem:exponential-polynomial}. In the following Lemma, we prove that such a polynomial exists.
\begin{restatable}[Polynomial approximation to $e^{-tx^2}$]{lemma}{lempolyapproxgaussian}
\label{lem:polynomial-gaussian-operator}
Suppose $x\in [-1,1]$, $\varepsilon\in [0,1/2)$ and $t\in \mathbb{R}^+$. Furthermore, suppose $d=\lceil\max\{te^2/2,\ln(2/\varepsilon)\}\rceil$. Then, there exists a polynomial $\tilde{q}_{t,d,d'}(x)$ of degree
d'=\lceil\sqrt{2d\ln(4/\varepsilon)}\rceil\in O\left(\sqrt{t}\log(1/\varepsilon)\right),
for which the following holds
\sup_{x\in [-1,1]}\abs{e^{-tx^2}-\tilde{q}_{t,d,d'}(x)}\leq \varepsilon.
\end{restatable}
\begin{proof}
This is proven in Sec.~\ref{sec:appendix-poly-approx} of the Appendix
\end{proof}

We use lemma \ref{lem:polynomial-gaussian-operator}, to obtain a block encoding of $e^{-tH^2}$, given an approximate block encoding of $H$. Subsequently, we shall show that this results in a robust quantum algorithm for preparing the ground state of $H$ under the assumptions we have considered.

\begin{lemma}
\label{lem:implementing-gaussian-operator}
Let $H$ be a Hermitian matrix with eigenvalues in $[-1,1]$ and $\varepsilon\in (0,1/2)$. Furthermore, suppose $t\in \mathbb{R}^{+}$ and $U_H$ is an $(1,a,\delta)$-block encoding of $H$, implementable in time $T_H$. Also, let $d=\lceil\max\{te^2/2,\ln(4/\varepsilon)\}\rceil$ and $d'=\lceil\sqrt{2d\ln(8/\varepsilon)}\rceil$. Then, provided
\delta\leq \dfrac{\varepsilon^2}{128 d~\ln(8/\varepsilon)},
%\delta\leq \dfrac{\varepsilon}{4\sqrt{d\log(8/\varepsilon)}}.
we can  implement an $(1,a+1,\varepsilon)$-block encoding of $e^{-t H^2}$ in cost
T=O\left(T_H \sqrt{t} \log(1/\varepsilon)\right).
\end{lemma}
\begin{proof}
Now suppose $\tilde{H}=\left(\bra{0}^{\otimes a}\otimes I\right)U_H \left(\ket{0}^{\otimes a}\otimes I\right)$. Then, from the definition of block encoding of operators, $\norm{H-\tilde{H}}\leq \delta$. Also, from Lemma \ref{lem:polynomial-gaussian-operator}, we can use the polynomial of degree $d'=\lceil\sqrt{2d\log(8/\varepsilon)}\rceil$ to implement $(1,a+1,\varepsilon/2)$-block encoding of $\tilde{q}_{t,d,d'}(\tilde{H})$ in cost
T=d'\in O\left( T_H \sqrt{t} \log(1/\varepsilon)\right).
The number of ancilla qubits increased by one because of the QSVT procedure. So, we have
\begin{align}
\norm{e^{-tH^2}-\tilde{q}_{d,d'}\left(\tilde{H}\right)}
&\leq \norm{e^{-tH^2}-\tilde{q}_{d,d'}\left(\tilde{H}\right)}+\norm{\tilde{q}_{d,d'}\left(\tilde{H}\right)-\tilde{q}_{d,d'}\left(H\right)}\\
&\leq \varepsilon/2+4d'\sqrt{\delta}\\
&\leq \varepsilon/2+\varepsilon/2=\varepsilon.~~~~~~~~~~[\text{~Substituting the value of $\delta$ and $d'$}].
\end{align}
\end{proof}
Now that we have a procedure to implement a block encoding of $e^{-tH^2}$, given an approximate block encoding of $H$, we can use this to obtain a circuit model quantum algorithm for preparing the $0$-eigenstate of $H$. As before, let us make some assumptions on the spectrum of $H$. We assume that we are given a Hamiltonian $H$ of unit norm with ground energy, $\lambda_0$ and we intend to prepare a state that is close to its ground state, $\ket{v_0}$. We assume that the gap between the ground state and the rest of the spectrum is lower bounded by $\Delta$. We also assume that we have knowledge of $E_0$ such that
|\lambda_0-E_0|\leq O\left(\Delta/\sqrt{\log\frac{1}{\eta\varepsilon}}\right).
\begin{lemma}
\label{lem:gsp-discrete}
Let $\varepsilon\in (0,1/2)$ and $H$ be a Hamiltonian. Furthermore, suppose we are given $U_H$, which is a $(1,a,\delta)$-block encoding of $H$, implemented in time $T_H$. Let $\ket{v_0}$ be the ground state of $H$ with eigenvalue $\lambda_0$ such that the value of $\lambda_0$ is known up to precision
$\eps_g\in \Oo\left(\Delta/\sqrt{\log\frac{1}{\eta\eps}}\right)$, where $\Delta$ is a lower bound on the spectral gap of $H$.

Additionally, let us assume access to a state preparation procedure $B$ which prepares a state $\ket{\psi_0}$ in time $T_{\psi_0}$ such that $|\braket{\psi_0|v_0}|\geq \eta$. Also, let
\delta\leq \dfrac{\varepsilon^2\eta^2}{512 d~\ln\left(\dfrac{16}{\eta\varepsilon}\right)},
%\delta\leq \dfrac{\varepsilon\eta}{8\sqrt{2d\log\left(\dfrac{16}{\eta\varepsilon}\right)}},
where, $d=\lceil\max\{te^2/2, \ln(8/\varepsilon)\}\rceil$, and
Then there exists a quantum algorithm that prepares a quantum state that is $O(\varepsilon)$-close to $\ket{v_0}$ with $\Omega(1)$ probability in cost
\begin{equation}
\label{eq:complexity-gsp}
\end{equation}
\end{lemma}
\begin{proof}
In lemma \ref{lem:implementing-gaussian-operator}, we replace $\varepsilon$ with $\varepsilon\eta/2$ to prepare an $(1,a+1,\varepsilon\eta/2)$-block encoding of $e^{-tH^2}$. Furthermore, we choose
\begin{equation}
\label{eq:value-of-t}
t\geq \dfrac{1}{2\Delta^2}\log\left(\dfrac{4(1-\eta^2)}{\eta^4\varepsilon^2}\right)=O\left(\dfrac{1}{\Delta^2}\log\left(\dfrac{1}{\eta\varepsilon}\right)\right).
\end{equation}
To get an $\varepsilon\eta/2$-precision in the block encoding the degree of the polynomial $\tilde{q}_{t,d,d'}(H')$ is
where $d=\left\lceil\max\{te^2/2, \ln\left(\frac{8}{\varepsilon\eta}\right)\}\right\rceil$.
This yields that the precision of block encoding of $H$ needs to be at least $\delta$-precise where,
\delta\leq \dfrac{\varepsilon^2\eta^2}{512d\ln\left(\dfrac{16}{\eta\varepsilon}\right)}.
Thus, with cost,
we prepare a quantum state that is $O\left(\varepsilon\eta/2\right)$-close to
\ket{\eta_t}=\ket{\bar{0}}\dfrac{e^{-tH^2}}{\sqrt{\braket{\psi_0|e^{-2tH^2}|\psi_0}}}\ket{\psi_0}+\ket{\Phi^\perp}.
Post-selected on obtaining $\ket{\bar{0}}$ in the first register, we obtain a quantum state that is $O(\varepsilon\eta/2)$-close to
\begin{equation}
\label{eq:state-phi}
\ket{\phi}=\dfrac{e^{-tH^2}}{\sqrt{\braket{\psi_0|e^{-2tH^2}|\psi_0}}}\ket{\psi_0},
\end{equation}
with amplitude $\sqrt{\braket{\psi_0|e^{-2tH^2}|\psi_0}}=\Omega(\eta)$, where the lower bound is obtained from Eq.~\eqref{eq:succ-prob-gsp-2}. Now by choosing $t$ as in Eq.~\eqref{eq:value-of-t} and replacing it in Eq.~\eqref{eq:distance-ground-state}, we have
\norm{\ket{v_0}-\ket{\phi}}\leq O(\varepsilon\eta/2).
By the triangle inequality, this implies, that the quantum state prepared is $O(\varepsilon\eta)$-close to $\ket{v_0}$ with probability $\Omega(\eta)$. So by using $O(1/\eta)$-rounds of amplitude amplification, we obtain a quantum state that is $O(\varepsilon)$-close to $\ket{v_0}$ with probability $\Omega(1)$. The overall cost will be
\end{proof}

Now suppose we have a generic Hamiltonian with eigenvalues in $[-1,1]$ such that the ground energy is unknown. Then the function $e^{-tH^2}$ helps prepare the $0$-eigenstate of $H$, even if this is not necessarily the ground state of $H$. One application of this procedure is that it leads to an optimal quantum linear systems algorithm as we shall see in the next section.
\\~\\
\section{Applications to Quantum linear systems}
\label{sec:qls}
The quantum linear systems algorithm can be stated as follows: Given access to a Hermitian matrix $H$ and some initial state $\ket{b}$, prepare the quantum state $\ket{x}=H^{-1}\ket{b}/\norm{H^{-1}\ket{b}}$. Ever since the first quantum algorithm for this problem by Harrow, Hassidim, and Lloyd [62], the quantum linear system algorithm has been widely studied. The complexity of this algorithm has been progressively improved through a series of results [25, 26, 46]. Recently, adiabatic-inspired approaches have also been reported [52, 53, 55]. Just like in the previous section, we apply ``Analog LCU'' to develop two quantum linear systems algorithms in continuous-time (Sec \ref{subsec:analog-qls}): the first one is an analog variant of the direct approach in [25] while the second one is more amenable to near-term implementations. Following this, we use the ``Single-Ancilla'' approach to develop two randomized quantum algorithms for this problem that are implementable on early fault-tolerant quantum computers (Sec.~\ref{subsec:single-ancilla-qls}).  Finally, just like in the previous section, we provide an algorithm for solving quantum linear systems using QSVT on fully fault-tolerant quantum computers (Sec.~\ref{subsec:qls-qsvt}).

Let us begin by formally stating the quantum linear systems problem.
\textbf{Quantum linear systems:~} Suppose we have access to a Hermitian matrix $H\in\mathbb{R}^{N\times N}$ such that its eigenvalues lie in the interval $[-1,-1/\kappa]\cup[1/\kappa, 1]$. Then, given a procedure that prepares the $N$-dimensional quantum state $\ket{b}$, a quantum linear systems algorithm prepares a quantum state that is $O(\varepsilon)$ - close to
\ket{x}=\dfrac{H^{-1}\ket{b}}{\norm{H^{-1}\ket{b}}}.
It is worth noting that the quantum linear systems algorithm is different from its classical counterpart in that by preparing $\ket{x}$, one does not have access to the entries of the classical vector $\vec{x}$. In quantum linear systems, thus, often one is interested in extracting useful information out of the state $\ket{x}$, such as estimating the expectation value $\braket{x|O|x}$, for some observable $O$. The assumption that $H$ is a Hermitian matrix is without loss of generality. Given any non-Hermitian $H\in \mathbb{R}^{N\times d}$, there exist efficient procedures to obtain a Hermitian matrix $\tilde{H}$ of dimension $(N+d)\times (N+d)$, such that the absolute value of its (non-zero) eigenvalues are equal to the non-zero singular values of $H$. Then, one may instead implement quantum linear systems with $\tilde{H}$ instead of $H$.

\subsection{Applying Analog LCU: Continuous-time quantum linear systems algorithms}
\label{subsec:analog-qls}
In this section, we develop analog quantum algorithms for solving quantum linear systems. Following the exposition in Sec.~\ref{subsec:analog-LCU}, we shall assume that we are given a system Hamiltonian $H$. We couple this Hamiltonian (the primary system) to two ancillary continuous-variable systems via the interaction Hamiltonian
\begin{equation}
\label{eq:interaction-hamiltonian}
H'=H\otimes \hat{y} \otimes \hat{z}.
\end{equation}
The primary system will be initialized in the quantum state $\ket{b}$ while the two ancillary systems will be in some continuous-variable states. The quantum algorithms developed in this subsection involve evolving the overall system according to $H'$ for some time. Following this, we shall show that the primary system is in the state $\ket{x}$ (or close to it) with an amplitude of $\Omega(1/\kappa)$.

We begin with the first quantum algorithm, which is an analog analogue of the quantum linear systems algorithm of [25].
\textbf{Continuous-time quantum linear systems algorithm:~} Consider the function $f(y)=ye^{-y^2/2}$, where $y\in\mathbb{R}$. As
\begin{align}
\implies & \int_{0}^{\infty}dy~f(xy)=1/x,
\end{align}
which holds for any $x\neq 0$. For any function $g(y)$, suppose its Fourier transform is $\mathcal{F}(g(y))=F(\omega)$, then $\mathcal{F}(g'(y))=i\omega F(\omega)$. If $g(y)=e^{-y^2/2}$, we have that $g'(y)=-ye^{-y^2/2}=-f(y)$. This implies,
\begin{align}
\label{eq:odd-function-fourier}
\dfrac{i}{\sqrt{2\pi}}\int_{-\infty}^{\infty} dz~ ze^{-z^2/2} e^{-izy}= ye^{-y^2/2},
\end{align}
\dfrac{1}{x}=\dfrac{i}{\sqrt{2\pi}}\int_{0}^{\infty} dt \int_{-\infty}^{\infty}dz~ ze^{-z^2/2} e^{-izxt}.
Next, we will prove via a lemma that the upper limit of the outer integral can be truncated at $T=\widetilde{O}(\kappa)$, without introducing significant error.
\begin{lemma}
\label{lem:truncation-2}
Suppose $\varepsilon >0, z\in \mathbb{R}$, and $x\in \mathbb{R}\setminus\{0\}$. Then there exists $T\in\Theta\left(\kappa\sqrt{\log(\kappa/\varepsilon)}\right)$, such that on the domain  $[-1,-1/\kappa]\cup[1/\kappa,1]$,
\begin{equation}
\label{eq:error-bound}
\left|\dfrac{1}{x}-\dfrac{1}{\sqrt{2\pi}}\int_{0}^{T} dt~\int_{-\infty}^{\infty}dz~ ze^{-z^2/2} e^{-izxt}\right|\leq \varepsilon.
\end{equation}
\end{lemma}
\begin{proof}
We have to evaluate the quantity
\left|\dfrac{1}{\sqrt{2\pi}}\int_{T}^{\infty}dt~\int_{-\infty}^{\infty}dz~ ze^{-z^2/2} e^{-izxt}\right|
We first evaluate the outer integral and obtain,
\begin{align}
\left|\dfrac{1}{\sqrt{2\pi}}\int_{T}^{\infty}dt~\int_{-\infty}^{\infty}dz~ ze^{-z^2/2} e^{-izxt}\right|&=\left|\int_{T}^{\infty}dt~ xt~e^{-x^2t^2/2}\right|~~~~~~[\text{~Using Eq.~\eqref{eq:odd-function-fourier}}~]\\
&=\left|\dfrac{1}{x}\int_{x^2T^2/2}^{\infty} dy~e^{-y}\right|~~~~~~~~~~~[~y=x^2t^2/2~]\\
&=\left|\dfrac{1}{x}\cdot e^{-x^2T^2/2}\right|\\
&\leq \dfrac{1}{|x|}\left|e^{-x^2T^2/2}\right|.
\label{eq:bound-gaussian}
\end{align}
Now for $T=\kappa \sqrt{2\log(\kappa/\varepsilon)}$, we have $\left|e^{-x^2T^2/2}\right|\leq \varepsilon/\kappa$. Now as $|x|\geq 1/\kappa$, we have that Eq.~\eqref{eq:bound-gaussian} is upper bounded by $\varepsilon$. So finally,
\begin{equation}
\left|\dfrac{1}{x}-\dfrac{1}{\sqrt{2\pi}}\int_{0}^{T} dt~\int_{-\infty}^{\infty}dz~ ze^{-z^2/2} e^{-izxt}\right|=\left|\dfrac{1}{\sqrt{2\pi}}\int_{T}^{\infty}dt~\int_{-\infty}^{\infty}dz~ ze^{-z^2/2} e^{-izxt}\right|\leq \varepsilon.
\end{equation}
\end{proof}

In order to design the analog quantum algorithm, consider that the effective interaction Hamiltonian is $H'=H\otimes \hat{y}\otimes \hat{z}$. While the system Hamiltonian $H$ is prepared in some input state $\ket{b}$, the first ancilla system is prepared in the first-excited state of a one-dimensional quantum Harmonic oscillator
\begin{equation}
\label{eq:first-excited-state-harmonic}
\ket{\psi_h}=\dfrac{1}{(2\pi)^{1/4}}\int_{-\infty}^{\infty} dy~ye^{-y^2/4}\ket{y}.
\end{equation}
The second ancilla system is in the ground state of a ``particle in a ring'' of diameter $1$, given by
\begin{equation}
\ket{\tau}=\int_{0}^{1}dz~\ket{z}.
\end{equation}
Then evolving the overall system according to $H'$ for time $T$, we obtain
\begin{align}
\ket{\eta_t}&=e^{-i\tilde{H}T}\ket{b}\ket{\psi_h}\ket{\tau}\\
&=\int_{0}^{1}dz~\int_{-\infty}^{\infty} \dfrac{dy}{(2\pi)^{1/4}}~y e^{-y^2/4} e^{-iyzHT}\ket{b}\ket{y}\ket{z}\\
&=\dfrac{1}{T}\int_{0}^{T}dt~\int_{-\infty}^{\infty} \dfrac{dz}{\sqrt{2\pi}}~ze^{-z^2/2} e^{-iztH}\ket{b}\ket{\psi_h}\ket{\tau}+\ket{\Phi}^{\perp}~~~~[~\text{Change of variable $t=Ty$}~]
\end{align}
Now, by choosing time $T=\Theta\left(\kappa\sqrt{\log(\kappa/\varepsilon)}\right)$, from Lemma \ref{lem:truncation-2}, we obtain a quantum state that is $O(\varepsilon/T)$-close to
\begin{align}
\ket{\eta_t}=\dfrac{H^{-1}}{T} \ket{b}\ket{\psi_h}\ket{\tau}+\ket{\Phi}^{\perp}.
\end{align}
By post-selecting on $\ket{\psi_h}$ in the second register, one obtains $\ket{x}$ in the first register with probability $\widetilde{\Omega}(1/\kappa^2)$. Alternatively, $\widetilde{O}(\kappa)$-rounds of amplitude amplification (which is a circuit model procedure) can yield a quantum state $O(\varepsilon)$-close to $\ket{x}$.

Although this procedure works in general, the quantum state $\ket{\tau}$ might be difficult to prepare experimentally. In fact, for continuous-variable systems, Gaussian states are the easiest to prepare and manipulate [49]. So, next, we provide a quantum algorithm for which it suffices to prepare both the ancillary registers in Gaussian states.
\textbf{Continuous-time quantum linear systems algorithm using only Gaussian states:~} The previous quantum algorithm sections require to prepare the non-Gaussian continuous-variable state
\begin{equation*}
\ket{\tau}=\dfrac{1}{\sqrt{T}}\int_{0}^{T}dz~\ket{z}.
\end{equation*}
Since Gaussian states are typically easier to generate and manipulate, let us design alternative algorithms using Gaussian states only. The general idea is to approximate $\int_{-\infty}^{+\infty}dt$ by $\int_{-\infty}^{+\infty}dt~e^{-t^2/2T^2}$ (rather than $\int_{-T}^{T}dt$) for large enough $T$. The analogue of Lemma~\ref{lem:truncation-2} becomes
\begin{lemma}
\label{lem:truncation-gaussian}
Suppose $\varepsilon >0, z\in \mathbb{R}$, and $x\in \mathbb{R}\setminus\{0\}$. Then, there exists $T\geq\kappa^{3/2}/\sqrt{\varepsilon}$, such that on the domain  $[1/\kappa,1]$,
\begin{equation}
\label{eq:truncated-inverse-gaussian}
\left|\dfrac{1}{x}-\dfrac{1}{2\pi}\int_{-\infty}^{+\infty} dt~e^{-t^2/2T^2}~\int_{-\infty}^{+\infty} dz~ e^{-z^2/2}e^{-ixtz} \right|\leq \Theta(\varepsilon).
\end{equation}
\end{lemma}

\begin{proof}
We have, using the fact that the Fourier transform of a Gaussian is a Gaussian
\begin{align*}
\dfrac{1}{2\pi}\int_{-\infty}^{+\infty} dt~e^{-t^2/2T^2}~\int_{-\infty}^{+\infty} dz~ e^{-z^2/2}e^{-ixtz}
&=\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty} dt~e^{-t^2/2T^2}~ e^{-x^2t^2/2}\\
&=\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty} dt~e^{-\left(x^2+1/T^2\right)t^2/2}
=\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty} dt~e^{-\tilde{x}^2t^2/2}\\
\end{align*}
where, we have set $\tilde{x}=\sqrt{x^2+1/T^2}$. Therefore, it remains to bound
\begin{align*}
\left|\dfrac{1}{x}-\dfrac{1}{\tilde{x}}\right|=\left|\dfrac{1}{x}\left(1-\dfrac{x}{\tilde{x}}\right)\right|
\leq\dfrac{1}{\abs{x}}\left|1-\dfrac{1}{\sqrt{1+1/x^2T^2}}\right|\leq \dfrac{1}{\abs{x}}\cdot\dfrac{1}{x^2T^2}\leq\varepsilon
\end{align*}
\end{proof}
Unfortunately, the scaling of $T$ is worse than for the non-Gaussian approach since $T$ scales as $\kappa^{3/2}$ (instead of linear)  and the dependence of precision is $1/\sqrt{\varepsilon}$ (rather than inverse-logarithmic). Moreover, as the Gaussian function is even, the procedure works only for positive semi-definite Hamiltonians. Nevertheless, this allows us to design a quantum linear systems algorithm using only Gaussian states as ancillae.

Let us again consider the interaction Hamiltonian $H'=H\otimes \hat{y}\otimes \hat{z}$, where $H$ is now some positive semidefinite Hamiltonian with its eigenvalues lying in $[1/\kappa,1]$. We prepare both the ancilla registers in a Gaussian state which, similarly to Sec.~\ref{subsec:analog-gsp}, is defined as follows
\begin{equation*}
\ket{\psi_g}=\dfrac{1}{(2\pi)^{1/4}} \int_{-\infty}^{\infty} dz~ e^{-z^2/4} \ket{z}.
\end{equation*}
Indeed, it suffices to let the state $\ket{b}\ket{\psi_g}\ket{\psi_g}$ evolve under Hamiltonian $H'$ for time $T$ to obtain
\begin{align*}
e^{-iH'T}\ket{b}\ket{\psi_g}\ket{\psi_g}=\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}dz~\int_{-\infty}^{\infty} dy~e^{-(y^2+z^2)/4} e^{-iyzHT}\ket{b}\ket{y}\ket{z}.
\end{align*}

If we choose some $T\geq\kappa^{3/2}/\sqrt{\varepsilon}$, we have
\begin{align}
\left(I\otimes\ket{\psi_g}\bra{\psi_g}\otimes \ket{\psi_g}\bra{\psi_g}\right)e^{-i\tilde{H}T}\ket{b}\ket{\psi_g}\ket{\psi_g}
&=\dfrac{1}{2\pi}\int_{-\infty}^{+\infty}dz~\int_{-\infty}^{+\infty} dy~e^{-(y^2+z^2)/2} e^{-iyzHT}\ket{b}\ket{\psi_g}\ket{\psi_g}\\
&=\dfrac{1}{2\pi T}\int_{-\infty}^{+\infty}dt~e^{-t^2/2T^2}~\int_{-\infty}^{+\infty} dz~e^{-z^2/2} e^{-itzH}\ket{b}\ket{\psi_g}\ket{\psi_g}
\end{align}
where we have used the change of variable $t=Ty$.
\begin{align}
e^{-i\tilde{H}T}\ket{b}\ket{\psi_g}\ket{\psi_g}&=\dfrac{1}{2\pi T}\int_{-\infty}^{+\infty}dt~e^{-t^2/2T^2}~\int_{-\infty}^{+\infty} dy~e^{-z^2/2} e^{-itzH}\ket{b}\ket{\psi_g}\ket{\psi_g}+\ket{\phi}^{\perp}\\
&=\dfrac{H^{-1}}{T} \ket{b}\ket{\psi_g}\ket{\psi_g}+\ket{\Phi}^{\perp}+O(\varepsilon/T)~~~~~~ [~\text{From Lemma \ref{lem:truncation-gaussian}}~].
\end{align}
So, ~by post-selecting on obtaining $\ket{\psi_g}$in the second  register, we obtain a state that is $O\left(\varepsilon/\kappa\right)^{3/2}$ -close to
\ket{x}=\dfrac{H^{-1}\ket{b}}{\norm{H^{-1}\ket{b}}},
with amplitude $\tilde{\Omega}\left(\sqrt{\varepsilon}/\kappa^{3/2}\right)$. Although the complexity of this algorithm is worse than the continuous-time quantum algorithm in the previous section, it requires only Gaussian states. Consequently, it is more suitable for being implementable in the near term.

One can improve the complexity of this quantum linear systems algorithm by replacing the Gaussian state in the second register with the flat state $\ket{\tau}$. For positive semidefinite Hamiltonians, if the first ancillary system is in a Gaussian state while the second one is in $\ket{\tau}$, we can still obtain a quantum state that is $O(\varepsilon/\kappa)$-close to the solution of the quantum linear systems in time $\widetilde{O}(\kappa)$. This follows from observing
\begin{equation}
\dfrac{1}{x}=\int_{-\infty}^{\infty} \dfrac{dt}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\dfrac{dz}{\sqrt{2\pi}} e^{-z^2/2}e^{-ixtz}=2\int_{0}^{\infty} \dfrac{dt}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\dfrac{dz}{\sqrt{2\pi}} e^{-z^2/2}e^{-ixtz}.
\end{equation}
We can truncate the outer integral to $T=\Theta(\kappa\sqrt{\log(\kappa/\varepsilon)})$, and introduce only a $\varepsilon$ error.

Finally, we believe several existing quantum technological platforms might already be able to engineer these interactions for system Hamiltonians of small dimensions. Our analog approach provides a more physically motivated model for implementing quantum linear systems. Next, we move on to the problem of sampling from the solution of quantum linear systems on early fault-tolerant quantum computers. For this, we make use of the ``Single-Ancilla LCU'' technique.
\subsection{Applying Single-Ancilla LCU: Sampling from the solution of quantum linear systems}
\label{subsec:single-ancilla-qls}

In the ``Single-Ancilla'' framework, we provide two randomized quantum algorithms for sampling from the solution of quantum linear systems. For the first algorithm, we consider the discrete LCU decomposition of $H^{-1}$ of Ref.~[25]. For the second method, we will draw inspiration from the adiabatic approaches to solve this problem [52, 53, 54, 55].
\textbf{Method 1: Direct implementation of $H^{-1}$}
\\~\\
In this approach, we will use the discrete version of the LCU decomposition of $H^{-1}$ by Childs, Kothari, and Somma [25]. We will directly use their results here. We begin by stating the discretized version of the expression in Lemma \ref{lem:truncation-2}.
\begin{equation}
\label{eq:discrete-lcu-inverse}
g(x)=\dfrac{i}{\sqrt{2\pi}}\sum_{j=0}^{J-1}\Delta_y\sum_{k=-K}^{K}\Delta_z z_k e^{-z^2_k/2}e^{-xy_jz_k},
\end{equation}
where $y_j=j\Delta_y$ and $z_k=k\Delta_z$, for some $J\in \Theta(\frac{\kappa}{\varepsilon}\log(\kappa/\varepsilon))$, $\Delta_y=\Theta(\varepsilon/\sqrt{\log(\kappa/\varepsilon)})$ and $\Delta_z=\Theta((\kappa\sqrt{\log(\kappa/\varepsilon)})^{-1})$. Then, Childs et al.~[25] proved that $\left|1/x-h(x)\right|\leq \varepsilon$ in the domain $[-1,-1/\kappa]\cup [1/\kappa,1]$. From this LCU it is clear that in order to approximate $H^{-1}$ in this domain, the time parameter of the Hamiltonian simulation is at most,
\begin{equation}
\label{eq:max-time-ham-sim-inversion}
t=\Theta(y_J z_K)=\Theta\left(\kappa\log(\kappa/\varepsilon)\right).
\end{equation}
Furthermore, from [25], the $\ell_1$-norm of the LCU coefficients were shown to be upper bounded by $\norm{c}_1=\Theta(\kappa\sqrt{\log(\kappa/\varepsilon)})$.

The normalization factor $\ell=\norm{H^{-1}\ket{b}}$ is important as we have seen in Sec.~\ref{subsec:gsp-single-ancilla}. In this case $1\leq \ell \leq \kappa$. As before, given an approximate estimate $\tilde{\ell}$ of $\ell^2$, for any observable $O$, we would want Algorithm \ref{algo:randomized-time-evolution} to output some $\mu$ such that
\begin{equation}
\left|\dfrac{\mu}{\tilde{\ell}}-\braket{x|O|x}\right|\leq \varepsilon.
\end{equation}
We show via the following Theorem that this is indeed the case for appropriate choices of $T$ and $\tau_{\max}$.
\begin{theorem}[Sampling from the solution of quantum linear systems]
\label{thm:inversion-estimation-2}
Let $H$ be a Hermitian matrix such that its non zero eigenvalues lie in $[-1,-1/\kappa]\cup [1/\kappa,1]$. Suppose there exists a procedure that prepares the quantum state $\ket{b}$ in cost $\tau_b$ and let $O$ be an observable. Define $\ell=\norm{H^{-1}\ket{b}}$ and let $\varepsilon, \delta, \gamma \in (0,1)$. Furthermore, assume that we know $\tilde{\ell}$, an approximation $\ell^2$ such that
\left|\tilde{\ell}-\ell^2\right|\leq \dfrac{\varepsilon \ell^2}{6\norm{O}}
suppose, $\varepsilon,\delta,\gamma\in (0,1)$ be parameters and $O$ be an observable. Then if
\gamma \leq \dfrac{\varepsilon \ell^2}{18\norm{O}\norm{H^{-1}}},
such that
\norm{H^{-1}-g(H)}\leq \gamma,
T\geq O\left(\dfrac{\norm{O}^2\kappa^4\log^2\left(\frac{\norm{O}\kappa}{\varepsilon}\right)\ln(2/\delta)}{\varepsilon^2}\right).
then Algorithm \ref{algo:randomized-time-evolution} allows outputs, with probability at least $1-\delta$, a parameter $\mu$ such that
\left|\dfrac{\mu}{\tilde{\ell}}-\braket{x|O|x}\right|\leq \varepsilon,
using $T$ calls to a Hamiltonian simulation procedure, where the cost of each such repetition is at most
\tau_{max}=O\left(\kappa\log\left(\dfrac{\norm{O}\kappa}{\varepsilon}\right)\right).
and only one ancilla qubit.
\end{theorem}
\begin{proof}
First observe that we have
\norm{g(H)-H^{-1}}\leq \gamma,
\gamma \leq \dfrac{\varepsilon}{18\kappa^2 \norm{O}},
which implies in fact,
\gamma \leq \dfrac{\varepsilon \ell^2}{18\norm{O}\norm{H^{-1}}},
as $\ell\geq 1$ and $\norm{H^{-1}}\leq \kappa$.
If $\rho_b=\ket{b}\bra{b}$, this upper bound in $\gamma$ ensures that from Theorem \ref{thm:randomized-time-evolution},
\left|\dfrac{\mu}{\ell^2}-\Tr[O~g(H)\rho_b g(H)^{\dag}]\right|\leq \varepsilon/6.
From Theorem \ref{thmsup:norm-approx}, and from the upper bound in $|\tilde{\ell}-\ell^2|$, this
\begin{equation}
\label{eq:bound-mu-by-l}
\left|\dfrac{\mu}{\tilde{\ell}}-\dfrac{\Tr[O~g(H)\rho_b g(H)^{\dag}]}{\ell^2}\right|\leq \varepsilon/2.
\end{equation}

Overall, Algorithm \ref{algo:randomized-time-evolution} outputs $\mu$ such that
\begin{align}
\left|\dfrac{\mu}{\tilde{\ell}}-\braket{x|O|x}\right|&\leq\left|\dfrac{\mu}{\ell^2}-\Tr[O~g(H)\rho_b g(H)^{\dag}]\right|+\left|\Tr[O~g(H)\rho_b g(H)^{\dag}-\braket{x|O|x}\right|
\label{eq:qls-distance-error}
\end{align}
The first expression in the RHS of Eq.~\eqref{eq:qls-distance-error} is bounded by $\varepsilon/2$ (from Eq.~\eqref{eq:bound-mu-by-l}). We bound the second expression to $\varepsilon/2$ from Theorem \ref{thm:randomized-time-evolution}. For this, the sample complexity is
\begin{equation}
T\geq O\left(\dfrac{\ln(2/\delta)\norm{O}^2\norm{c}^4_1}{\varepsilon^2 \ell^4}\right)=O\left(\dfrac{\norm{O}^2\kappa^4\log^2\left(\frac{\norm{O}\kappa}{\varepsilon}\right)\ln(2/\delta)}{\varepsilon^2}\right).
\end{equation}
The cost of each repetition is bounded by $\tau_{max}+\tau_b$, where $\tau_{\max}$ is determined by the upper bound in $\gamma$ as
\tau_{\max}= O\left(\kappa\log\left(\dfrac{\norm{O}\norm{H^{-1}}\kappa}{\varepsilon \ell^2}\right)\right)=O\left(\kappa\log\left(\dfrac{\norm{O}\kappa}{\varepsilon}\right)\right),
where we have used $\ell\geq 1$ and $\norm{H^{-1}}\leq \kappa$.
\end{proof}
For our next method, we draw inspiration from the problem of sampling from ground states. This approach has a slightly better sample complexity (up to log factors) but needs an extra ancilla qubit.
\textbf{Method 2: Quantum linear systems as a ground state preparation problem}
In Sec.~\ref{subsec:gsp-single-ancilla}, we have discussed that given a Hamiltonian $H$, the function $f(H)=e^{-tH^2}$ can be used to prepare the ground state of $H$ when $H$ is positive semi-definite or when the ground energy is known to a certain precision. However, for general Hamiltonians with eigenvalues in $[-1,1]$, $f(H)$ helps prepare the $0$-eigenstate of $H$. We shall exploit this fact in this section.

In fact, following the adiabatic inspired approaches for solving quantum linear systems in Refs~[52, 53, 54, 55], we know that, given any Hamiltonian $H$, one can construct a Hamiltonian $H'$ such that it's $0$-eigenstate encodes the solution to the linear systems problem. For simplicity, we restrict our attention to the case where $H$ is positive semi-definite, i.e.\ its eigenvalues are in $[1/\kappa,1]$ but this is without loss of generality. Our analysis holds for any $H$ as from any such $H$, one can construct $H'$. Furthermore, we consider that $H$ is accessible via a block encoding.

Consider the Hamiltonian
\begin{equation}
\label{eq:inversion-operator}
H'=\ket{0}\bra{1}\otimes H Q_b+\ket{1}\bra{0} \otimes Q_b H,
\end{equation}
where $Q_b=I-\ket{b}\bra{b}$. $H'$ can be constructed efficiently from an $(1, a, \varepsilon)$-block encoding of $H$. Let us consider the spectrum of $H'$ following Refs.~[52, 53, 54] and relate it to the spectrum of $H$. Once again, when $H$ is not positive definite, then there also exists a different encoding of $H'$ which can be applied (by adding two ancilla qubits). Moreover, the $0$-eigenstate of this Hamiltonian also encodes the solution to the quantum linear systems problem. For simplicity, we now restrict our attention to positive definite $H$.

It is easy to see that the null space of $H'$ is  $\mathrm{Span}\{\ket{0,x},\ket{1,b}\}$. In other words, the state $\ket{0,x}$ is a $0$-eigenstate of $H'$. Furthermore, if we start with the state $\ket{0,b}$, the dynamics never leaves the space spanned by states having $\ket{0}$ in the first register. Hence, $\ket{0,x}$ is the unique $0$-eigenstate of $H'$ in this subspace.

Also, one can verify that the minimum non-zero eigenvalue of $H'$ is at least $1/\kappa$ away from zero, i.e.
\Delta=\min_{\lambda\neq 0}|\lambda|=1/\kappa.

We will first show how to construct an $(1, a+2, \varepsilon)$-block encoding of $H'$ given a $(1, a, \varepsilon)$-block encoding of $H$, prepared with cost $T_H$. Moreover, if $T_b$ is the cost of implementing a procedure $B$, such that $B\ket{\bar{0}}=\ket{b}$, then the cost of preparing $H$ is $O(T_H+T_b)$. We formally prove this via the following lemma
\begin{lemma}
\label{lem:block encoding-ham}
Let $H$ be a Hermitian matrix such that its eigenvalues lie in $[1/\kappa,1]$. Suppose $U_H$, which is a $(1,a,\varepsilon)$-block encoding of $H$, can be implemented with cost $T_H$. Also, suppose we have access to the unitary $B$ such that $B\ket{\bar{0}}=\ket{b}$. Furthermore, define $Q_b=I-\ket{b}\bra{b}$. Then we can implement $U_{H'}$, which is an $(1,a+2,\varepsilon)$-block encoding of $H'=\ket{0}\bra{1}\otimes H Q_b+\ket{1}\bra{0}\otimes Q_b H$ in cost $O\left(T_H+T_b\right)$.
\end{lemma}
\begin{proof}
First, observe that $H'$ can be written as a product of three matrices. That is,
\begin{equation}
\label{eq:hamiltonian-breakdown}
I & 0\\
0 & Q_b
\end{pmatrix}
\begin{pmatrix}
0 & H\\
H & 0
\end{pmatrix}
\begin{pmatrix}
I & 0\\
0 & Q_b
\end{pmatrix}.
\end{equation}
Also, the matrix $Q_b$ can be written as
where $\Pi_0=\ket{0}\bra{0}$.

A $(1,1,0)$-block encoding of $Q_b$ is implemented by the circuit
V_Q= (I\otimes B)(H\otimes I)\left(I\otimes\ket{0}\bra{0}+e^{i\pi\Pi_0}\otimes \ket{1}\bra{1}\right)(H\otimes I) (I\otimes B^\dag).

The second matrix in the RHS of Eq.~\eqref{eq:hamiltonian-breakdown} is simply $H\otimes\sigma_x$. This can be implemented using $V_H=U_H\otimes \sigma_x$, which is a $(1,a,\varepsilon)$-block encoding of $H\otimes \sigma_x$. The first and the third matrix can be implemented by $\tilde{V}_Q=\ket{0}\bra{0}\otimes I +\ket{1}\bra{1}\otimes V_Q$, which is a $(1,2,0)$-block encoding of these matrices. Thus, we have the product of three block encodings. Then the unitary $U_H=(I\otimes \tilde{V}_Q)(I\otimes V_H)(I\otimes \tilde{V}_Q)$ is an $(1,a+2,\varepsilon)$ block encoding of $H'$, requiring cost $O(T_H+T_b)$.
\end{proof}

Note that $H$ need not be accessed via a block encoding. If $H$ is a linear combination of Pauli matrices, $H'$ can be built from $H$ using a short-depth quantum circuit. However, this requires one additional qubit which is an overhead. Thus, we have now mapped the quantum linear systems algorithm to prepare the zero-eigenstate of the Hermitian matrix $H'$. So in order to estimate $\braket{x|O|x}$, we can directly use Theorem \ref{thm:ground-state-sampling}, where the initial state is $\ket{0}\ket{b}$ and the $0$-eigenstate is $\ket{0}\ket{x}$.

Recall, using Algorithm \ref{algo:randomized-time-evolution} we can output $\mu$ such that for any observable $O$,
\begin{equation}
\left|\dfrac{\mu}{\ell^2} - \braket{x|O|x}\right|\leq \varepsilon,
\end{equation}
where $\ell=\norm{e^{-tH'^2}\ket{0,b}}\geq 1/\kappa$. However, we only need an approximation $\tilde{\ell}$ of $\ell^2$. We formally state our result via the following theorem
\begin{theorem}
\label{thm:inversion-estimation-1}
Let $H'$ be the Hamiltonian defined in Eq.~\eqref{eq:inversion-operator}, constructed from a given Hermitian matrix $H$ with its eigenvalues lying in $[1/\kappa,1]$. Suppose there exists a procedure that prepares the quantum state $\ket{b}$~in cost $\tau_b$. Let $O$ be some observable. Also for some,
define $\ell^2=\braket{0,b|e^{-2tH'^2}|0,b}$ . Suppose we know $\tilde{\ell}$, an approximation of $\ell^2$ such that
\left|\tilde{\ell}-\ell^2\right|\leq \dfrac{\varepsilon\ell^2}{6\norm{O}}.
Furthermore suppose, $\varepsilon,\delta,\gamma\in (0,1)$ be parameters. Then if
\gamma \leq \dfrac{\varepsilon}{36\kappa^2 \norm{O}},
T\geq O\left(\dfrac{\norm{O}^2\kappa^4\ln(2/\delta)}{\varepsilon^2}\right).
then Algorithm \ref{algo:randomized-time-evolution} allows outputs, with probability at least $1-\delta$, a parameter $\mu$ such that
\left|\dfrac{\mu}{\tilde{\ell}}-\braket{x|O|x}\right|\leq \varepsilon,
using $T$ calls to the quantum circuit in Fig.~\ref{fig:single-ancilla-circuit}, where the cost of each such repetition is at most $\tau_{\max}+\tau_b$, where
\tau_{max}=O\left(\kappa\log\left(\dfrac{\norm{O}\kappa}{\varepsilon}\right)\right).
\end{theorem}
\begin{proof}
First, observe that the upper bound in $\gamma$ ensures that
\gamma \leq \dfrac{\varepsilon \ell^2}{36\norm{O}}.
The rest of the proof follows from Theorem \ref{thm:ground-state-sampling} where we simply replace $\Delta=1/\kappa$ and $\eta\geq 1/\kappa$. So, we output such a $\mu/\tilde{\ell}$, with probability at least $1-\delta$, using
T\geq O\left(\dfrac{\norm{O}^2\kappa^4\ln(2/\delta)}{\varepsilon^2}\right)
samples of the quantum circuit in Fig.~\ref{fig:single-ancilla-circuit} where the cost of each such repetition is at most $\tau_{\max}+\tau_b$ where
\tau_{max}=O\left(\kappa\log\left(\dfrac{\norm{O}\kappa}{\varepsilon}\right)\right).
\end{proof}
Overall, the sample complexity of this method is slightly better than the randomized quantum algorithm in the first method. However, constructing $H'$ requires one extra qubit as compared to the previous approach. Nevertheless, this helps us connect between a ground state preparation procedure and the quantum linear systems algorithm which serves as an inspiration for our next quantum algorithm which uses QSVT to prepare the $0$-eigenstate of $e^{-tH'^2}$, given a block encoding of $H$.
\subsection{Quantum linear systems for fully fault-tolerant quantum computers}
\label{subsec:qls-qsvt}

The quantum algorithm is quite simple and exploits the connection observed in Method 2 of the previous section. Namely, the fact that given a block encoding of a Hermitian matrix $H$ whose eigenvalues are in $[-1, -1/\kappa]\cup [1/\kappa,1]$, each can efficiently construct a block encoding of a Hamiltonian $H'$ (See Lemma \ref{lem:block encoding-ham}) whose $0$-eigenstate is $\ket{0}\ket{x}$.

We use QSVT to implement a polynomial approximation of the Gaussian function (Lemma \ref{lem:implementing-gaussian-operator}), to the block encoding of $H'$. Formally,

\begin{lemma}
\label{lem:linear-systems-discrete}
Let $\varepsilon\in (0,1/2)$. Suppose $H$ is a Hermitian matrix with eigenvalues in $[1/\kappa,1]$ and,
\delta\in o\left(\dfrac{\varepsilon^2}{\kappa^4\log^2(\kappa/\varepsilon)}\right).
Furthermore, suppose that $U_H$ is an $(1,a,\delta)$-block encoding of $H$, implementable in cost $T_H$. Also, let $B$ be a unitary procedure that prepares $\ket{b}$ in time $T_b$. Then there exists a procedure that prepares a state that is $O(\varepsilon)$-close to,
\ket{x}=\dfrac{H^{-1}\ket{b}}{\norm{H^{-1}\ket{b}}}
with cost $O\left(\kappa^2 (T_H+T_b) \log(\kappa/\varepsilon) \right)$.
\end{lemma}
\begin{proof}
Given an $(1,a,\delta)$-block encoding of $H$, we can obtain, in time $O(T_H+T_b)$, an $(1, a+2,\delta)$-block encoding of the Hamiltonian $H'$ using Lemma \ref{lem:block encoding-ham}.

We now simply apply Lemma \ref{lem:gsp-discrete} by replacing $\Delta=1/\kappa$,
\eta=\braket{b|A^{-1}|b}\geq 1/\kappa.
and $T_H=O(T_A+T_b)$. Note that the precision in block encoding required is
\delta \in o\left(\dfrac{\varepsilon^2\eta^2}{t\log^2(1/\varepsilon\eta)}\right),
where $\sqrt{t}=O\left(\Delta^{-1}\log(\eta^{-1}\varepsilon^{-1})\right)$ and $\eta\geq 1/\kappa$. The complexity is obtained by the aforementioned substitutions in Eq.~\eqref{eq:complexity-gsp}.
\end{proof}
\section{Applications to quantum walks}
\label{sec:quantum-random-walks}

So far, we have seen applications of the ``Analog LCU'' and the ``Single-Ancilla LCU'' approaches. In this section, we will show that the ``Ancilla-free'' LCU can be applied to the framework of quantum walks.  Recall from Sec.~\ref{subsec:ancilla-free-LCU}, that this approach is useful when we are interested in the projection of the LCU state $f(H)\ket{\psi_0}$ in some subspace of interest. In such scenarios, it suffices to prepare an average density matrix $\rho$ by sampling the unitaries $U_j$ according to the distribution of the LCU coefficients. This is because the projection of $\rho$ in this subspace is at least as large (See Theorem \ref{thm:ancilla-free-LCU-theorem}). We will show this is precisely the case for spatial search by quantum walks.

We first discuss the optimal quantum spatial search algorithm by discrete-time quantum walks, for which we provide two quantum algorithms. The first relies on fast-forwarding discrete-time quantum walks [31]. For this algorithm, we formalize the unproven observation of Ref.~[41] - where the authors stated that the LCU could indeed be bypassed. Our second quantum algorithm relies on fast-forwarding continuous-time random walks, which also fits nicely in the ``Ancilla-free LCU framework''. For completeness, we also briefly outline the recent optimal spatial search algorithm by continuous-time quantum walk [32]. Finally, we show how one can obtain a discrete-time quantum walk from a continuous-time quantum walk (and vice versa) using the frameworks of block encoding and QSVT.