Quantum linear systems algorithms: a primer

TL;DR

This paper provides a detailed overview of the HHL quantum algorithm for solving linear systems, including its improvements, subroutines, and applications, highlighting its potential exponential speed-up over classical methods.

Contribution

It offers a comprehensive explanation of the HHL algorithm, its enhancements, and related quantum subroutines, serving as a foundational primer on quantum linear system solvers.

Findings

HHL achieves exponential speed-up over classical algorithms.

Improvements include variable-time amplitude amplification and Fourier/Chebyshev-based LCUs.

Discussion of quantum data loading and quantum singular value estimation.

Abstract

The Harrow-Hassidim-Lloyd (HHL) quantum algorithm for sampling from the solution of a linear system provides an exponential speed-up over its classical counterpart. The problem of solving a system of linear equations has a wide scope of applications, and thus HHL constitutes an important algorithmic primitive. In these notes, we present the HHL algorithm and its improved versions in detail, including explanations of the constituent sub- routines. More specifically, we discuss various quantum subroutines such as quantum phase estimation and amplitude amplification, as well as the important question of loading data into a quantum computer, via quantum RAM. The improvements to the original algorithm exploit variable-time amplitude amplification as well as a method for implementing linear combinations of unitary operations (LCUs) based on a decomposition of the operators using Fourier and Chebyshev series. Finally, we discuss a linear solver based on the quantum singular value estimation (QSVE) subroutine.

Figures

Tables

Table 1. Table 1: Runtime comparison between various quantum linear systems algorithms and the best general-purpose classical algorithm, conjugate gradient (CG). The parameter N 𝑁 N is the dimension of the system, s 𝑠 s is the sparsity of the matrix A 𝐴 A , κ 𝜅 \kappa is the condition number of A 𝐴 A , ϵ italic-ϵ \epsilon is the desired precision and ∥ ⋅ ∥ F subscript norm ⋅ 𝐹 \norm{\,\cdot\,}_{F} is the Frobenius norm.

Problem	Algorithm	Runtime Complexity
LSP	CG [She94]	$𝒪\left(Ns\kappa \log \left(1/ϵ\right)\right)$
QLSP	HHL [HHL09]	$𝒪\left(\log \left(N\right)s^2{\kappa }^2/ϵ\right)$
QLSP	VTAA-HHL [Amb10]	$𝒪\left(\log \left(N\right)s^2\kappa /ϵ\right)$
QLSP	Childs et. al. [CKS17]	$𝒪\left(s\kappa \mathrm{polylog}\left(s\kappa /ϵ\right)\right)$
QLSP	QLSA [WZP18]	$𝒪\left({\kappa }^2\mathrm{polylog}\left(n\right){\Vert A\Vert }_F/ϵ\right)$

Equations

$$|\psi ⟩=\sum _{i_1,\dots ,i_n}{\alpha }_{i_1\dots i_n}|i_1\dots i_n⟩,$$

$$p(m)=⟨\ast |\psi M_m^{†}{M}_m|\ast ⟩\psi$$

$$\frac{M_m|\ast ⟩\psi }{\sqrt{⟨\ast |\psi M_m^{†}{M}_m|\ast ⟩\psi }}.$$

$$|x⟩|q⟩\to |x⟩|q\oplus f(x)⟩$$

$$U_f=\sum _{x\in \{0,1{\}}^n}\sum _{q\in \{0,1{\}}^m}\text{outerproduct}\ast xx\otimes \text{outerproduct}\ast q\oplus f(x)q.$$

$$\text{QFT}:|\ast ⟩x\to \frac{1}{\sqrt{N}}\sum _{k=0}^{N-1}{\omega }^{x\cdot k}|\ast ⟩k.$$

$${\text{QFT}}^{†}:|\ast ⟩k\to \frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}{\omega }^{-k\cdot x}|\ast ⟩x.$$

$$|f_{j}\rangle=\frac{1}{\sqrt{N}}\sum_{k_{1},\ldots,k_{n}}\exp\big(i2\pi jk_{1}2^{-1}\big{missing})\exp\big(i2\pi jk_{2}2^{-2}\big{missing})\ldots\exp\big(i2\pi jk_{n}2^{-n}\big{missing})|k_{1}\ldots k_{n}\rangle.$$

$$|f_{j}\rangle=\frac{1}{\sqrt{N}}\Big{(}|0\rangle+\exp\big(i2\pi j2^{-1}\big{missing})|1\rangle\Big{)}\Big{(}|0\rangle+\exp\big(i2\pi j2^{-2}\big{missing})|1\rangle\Big{)}\ldots\Big{(}|0\rangle+\exp\big(i2\pi j2^{-n}\big{missing})|1\rangle\Big{)}.$$

$$|f_{j}\rangle=\frac{1}{\sqrt{N}}\Big{(}|0\rangle+\exp\big(i2\pi 0.j_{n}\big{missing})|1\rangle\Big{)}\Big{(}|0\rangle+\exp\big(i2\pi 0.j_{n-1}j_{n}\big{missing})|1\rangle\Big{)}\ldots\Big{(}|0\rangle+\exp\big(i2\pi j0.j_{1}\cdots j_{n}\big{missing})|1\rangle\Big{)}.$$

$$R_{k}=\left[{\begin{array}[]{cc}1&0\\ 0&e^{i2\pi\frac{1}{2^{k}}}\\ \end{array}}\right],$$

$$i\frac{d}{dt}|\ast ⟩\Psi =\hat{H}|\ast ⟩\Psi ,$$

$$e^{-iU\hat{H}{U}^{†}t}=U{e}^{-i\hat{H}t}{U}^{†},$$

$$|\ast ⟩i,0\to |\ast ⟩i,H_{ii}$$

$$\to e^{-i{H}_{ii}t}|\ast ⟩i,H_{ii}$$

$$\to e^{-i{H}_{ii}t}|\ast ⟩i,0=e^{-i\hat{H}t}|\ast ⟩i\otimes |\ast ⟩0.$$

$$\exp ({\hat{H}}_1+{\hat{H}}_2)$$

$$=$$

$$e^{-i({\hat{H}}_1+{\hat{H}}_2)t}=\underset{m\to \infty }{\lim }{\left(e^{-i{\hat{H}}_{1}t/m}{e}^{-i{\hat{H}}_{2}t/m}\right)}^m.$$

$$\text{norm}{e^{-i({\hat{H}}_1+{\hat{H}}_2)t}-{\left(e^{-i{\hat{H}}_{1}t/m}{e}^{-i{\hat{H}}_{2}t/m}\right)}^m}_2\le ϵ,$$

$${\left(e^{-i{\hat{H}}_{1}t/m}{e}^{-i{\hat{H}}_{2}t/m}\right)}^m$$

$$={\left(e^{-i({\hat{H}}_1+{\hat{H}}_2)t/m}+\mathcal{O}\left(\frac{t^2\max {({\hat{H}}_1,{\hat{H}}_2)}^2}{m^2}\right)\right)}^m.$$

$$(A+B{)}^m=\underset{0^{th}\text{ order in }B}{\underbrace{A^m}}+\underset{\text{m first order in }B}{\underbrace{A^{m-1}B+A^{m-2}BA+A^{m-3}B{A}^2++\dots +B{A}^{m-1}}}+$$

$$+\underset{\text{Second order in }B}{\underbrace{A^{m-2}{B}^2+A^{m-3}BAB+A^{m-3}{B}^{2}A}}+\dots +\underset{m^{th}\text{ order in }B}{\underbrace{B^m}},$$

$$e^{-i({\hat{H}}_1+{\hat{H}}_2)t(m-k-1)/m}\mathcal{O}\left(\frac{t^2\max {\left({\hat{H}}_1,{\hat{H}}_2\right)}^2}{m^2}\right)e^{-i({\hat{H}}_1+{\hat{H}}_2)tk/m},$$

$$e^{-i({\hat{H}}_1+{\hat{H}}_2)t}+\tilde{\mathcal{O}}\left(\frac{t^2\max {({\hat{H}}_1,{\hat{H}}_2)}^2}{m}\right).$$

$$\text{norm}{\left(e^{-i{\hat{H}}_{1}t/m}{e}^{-i{\hat{H}}_{2}t/m}\right)}^m-e^{-i({\hat{H}}_1+{\hat{H}}_2)t}=\mathcal{O}\left(\frac{t^2\max {\left(\text{norm}{{\hat{H}}_1}_2,\text{norm}{{\hat{H}}_2}_2\right)}^2}{m}\right).$$

$$e^{-i({\hat{H}}_1+\dots +{\hat{H}}_k)t}=\underset{m\to \infty }{\lim }{\left(e^{-i{\hat{H}}_{1}t/m}\cdots e^{-i{\hat{H}}_{k}t/m}\right)}^m.$$

$$\hat{B}=\left(\begin{array}[]{c c}0&\hat{H}\\ \hat{H}^{\dagger}&0\end{array}\right).$$

$$e^{-i({\sigma }_x\otimes \hat{H})t}|\ast ⟩+|\ast ⟩\psi =|\ast ⟩+\otimes e^{-i\hat{H}t}|\ast ⟩\psi .$$

$$|\ast ⟩a,b,h_a\to h_a|\ast ⟩a,b,h_a\to ,$$

Full text

1. 1 Introduction

2. 1.1 Motivation

3. 1.2 Quantum linear systems algorithms

4. 1.3 Quantum computing

5. 1.4 Quantum algorithms and machine learning

6. 1.5 Structure

7. 2 Quantum algorithms: fundamental components

8. 2.1 Notation and conventions

9. 2.2 Quantum Fourier transform

10. 2.2.1 The discrete Fourier transform

11. 2.2.2 The quantum Fourier transform

12. 2.2.3 Implementation of the QFT

13. 2.3 Hamiltonian simulation

14. 2.3.1 Trotter-Suzuki methods

15. 2.3.2 Graph-colouring method

16. 2.4 Erroneous Hamiltonian simulation

17. 2.4.1 Modern methods for Hamiltonian simulation

18. 2.5 Quantum phase estimation

19. 2.6 Phase kickback

20. 2.7 Amplitude amplification

21. 2.8 The uncompute trick

22. 2.9 Quantum RAM

23. 3 HHL

24. 3.1 Problem definition

25. 3.2 The HHL algorithm

26. 3.2.1 Algorithm summary

27. 3.2.2 Algorithm details

28. 3.3 Error analysis

29. 3.3.1 A short analysis

30. 3.3.2 Detailed analysis

31. 3.4 Matrix inversion is $\BQP$-complete

32. 3.5 Optimality of HHL

33. 3.6 Non-Hermitian Matrices

34. 4 Improvements to the HHL

35. 4.1 Variable time amplitude amplification for improved condition number dependence

36. 4.2 Exponentially improved precision number dependence

37. 4.3 QSVE based QLSA

38. 4.3.1 Singular values and subspaces

39. 4.3.2 The QSVE algorithm

40. 4.3.3 Extension to QLSA

1 Introduction

1.1 Motivation

Quantum computing was introduced in the $1980$s as a novel paradigm of computation, whereby information is encoded within a quantum system, as opposed to a system governed by the laws of classical physics. Wider interest in quantum computation has been motivated by Shor’s quantum algorithm for integer factorisation [Sho99], which provides an exponential speed-up over the best known classical algorithm for the same task. If implemented at scale, this would have severe security consequences for the ubiquitous RSA cryptographic protocol. Since then, a number of quantum algorithms demonstrating advantage over classical methods have been developed for a substantial variety of tasks; for a detailed survey, the reader is directed to  [Cle+98, Mon15].

Consistent advances on both theoretical and experimental research fronts have meant that the reality of a quantum computer has been edging ever closer. Quantum systems are extremely sensitive to noise, with a primary challenge being the development of error correction in order to achieve fault tolerance [ETC17]. Nonetheless, the current advances observed over recent years in the quest to build a scalable universal quantum computer from both academic and industrial groups raise the question of applications of such a device.

Although powerful quantum algorithms have been devised, their application is restricted to a few use cases. Indeed, the design of a quantum algorithm directly relies on exploiting the laws and features of quantum mechanics in order to achieve a speed-up. More precisely, in quantum computing, a quantum state is first prepared, to which quantum operations are applied before final measurements are performed, thus yielding classical data. As discussed in [Aar15], this model raises a number of challenges. In particular, careful consideration of how classical data can be input and obtained as output is crucial to maintaining the theoretical advantage afforded by quantum algorithms.

The question of solving a system of linear equations can be found at the heart of many problems with a wide scope of applications. An important result in recent years has been the Quantum Linear System algorithm (QLSA) [HHL09], also called Harrow-Hassidim-Lloyd (HHL) algorithm, which considers the quantum version of this problem. In particular, the HHL algorithm run on the quantum problem (that is, with quantum states as input and output) offers an exponential speed-up over the best known classical algorithm run on the classical problem. In the following, we present a complete review of the HHL algorithm and subsequent improvements for sampling from the solutions to linear systems. We have aimed to be as complete as possible, including relevant background where necessary. We assume knowledge of elementary linear algebra and some experience with analysis of classical algorithms.

1.2 Quantum linear systems algorithms

Solving a linear system is the task of taking a given matrix $A$ and vector $\mathbf{b}$ and returning a vector $\mathbf{x}$ satisfying $A\mathbf{x}=\mathbf{b}$. As a preview of what is coming up ahead, Table 1 shows the runtime of the best classical algorithm for solving linear systems, conjugate gradient (CG), compared with the quantum algorithms we shall introduce throughout these notes.

We note that CG solves a linear system completely, i.e. it returns the solution vector $\mathbf{x}$. The quantum algorithms allow one to sample from the solution efficiently, providing one has an efficient preparation method for the input state, i.e. a mapping from the vector $\mathbf{b}$ to a quantum state $\ket*{b}$.

We shall formally define the linear systems problem (LSP) and its quantum variant (QLSP) in section 3.1 and will discuss the differences between the two. We will discuss efficient state preparation in section 2.9.

1.3 Quantum computing

In this section, we introduce gate-model quantum computation, the computational model which will be used throughout. For a complete introduction to quantum computing we refer the reader to Nielsen and Chuang [NC02].

In classical computing, the input is a classical bit string which, through the application of a circuit, is transformed to an output bit string. This is achieved via the application of a finite number of classical gates picked from a universal gate set such as $\{$NAND$\}$. This framework is known as the classical circuit model of computation. In the case of quantum computing, there exist different yet equivalent frameworks in which the computation can be described: the quantum circuit model [NC02], measurement-based quantum computing (MBQC) [RBB03] and adiabatic quantum computing [Far+00]. In the following, we shall concentrate on the quantum circuit model, which is the closest quantum analogue of the classical circuit model as well as the model in which quantum algorithms are generally presented.

In classical computing the fundamental unit of information is the bit, which is either [math] or $1$, whereas in quantum computing, it is the qubit, $|\psi\rangle=\alpha_{0}|0\rangle+\alpha_{1}|1\rangle$, such that $\alpha_{0},\alpha_{1}\in\mathbb{C}$ and $|\alpha_{0}|^{2}+|\alpha_{1}^{2}|=1$. This is represented by a two-dimensional column vector belonging to a complex Hilbert space $\mathcal{H}\cong\mathbb{C}^{2}$. We shall denote by $\bra*{\psi}$ the conjugate-transpose of $\ket*{\psi}$. The states $|0\rangle$ and $|1\rangle$ are basis vectors corresponding to a bit value of ‘0’ and ‘1’ respectively, and we can write $\ket*{0}\cong\smallmatrixquantity(1\\ 0)$ and $\ket*{1}\cong\smallmatrixquantity(0\\ 1)$. Thus, we have that a qubit is a normalised complex superposition over these basis vectors. Multiple qubits are combined using the tensor product, that is, for two qubits $\ket*{\phi},\ket*{\psi}$, their joint state is given by $\ket*{\phi}\otimes\ket*{\psi}$. Thus, an $n$-qubit quantum state can be expressed as

[TABLE]

where $i_{k}\in\{0,1\}$, $\sum_{i_{1}\ldots i_{n}}|\alpha_{i_{1}\ldots i_{n}}|^{2}=1$ and $\ket*{i_{1}\ldots i_{n}}\equiv\ket*{i_{1}}\otimes\cdots\otimes\ket*{i_{n}}$. Note that $2^{n}$ complex coefficients are required to describe a quantum state, a number growing exponentially with the system’s size. We call the basis $\left\{\ket*{i_{1}\ldots i_{n}}\mid i_{k}\in\{0,1\}\right\}$ the computational basis, as each basis vector is described by a string of $n$ bits.

There are two types of operations we can perform on a quantum state: unitary operators and measurements. A unitary operator $U$ has the property that $UU^{\dagger}=U^{\dagger}U=I$, i.e. its inverse is given by the hermitian conjugate. Furthermore, this implies that the operator is norm-preserving, that is, unitary operators map quantum states to quantum states. A unitary operator on $n$ qubits can be expressed as matrix of dimension $2^{n}\times 2^{n}$. Moreover, we have that unitary operators are closed under composition. A measurement is described by a collection of (not necessarily unitary) operators $\{M_{k}\}$, where the index $k$ indicates a given measurement outcome. The operators $M_{k}$ act on the state’s Hilbert space and satisfy the completeness equation, $\sum_{k}M_{k}^{\dagger}M_{k}=I$. For a quantum state $\ket*{\psi}$, the probability of measuring outcome $m$ is given by

[TABLE]

and the resulting quantum state is then

[TABLE]

The completeness equation encodes the fact that measurement probabilities over all outcomes sum to unity. A computational basis measurement, $\{M_{x}\}$ for $x\in\{0,1\}^{n}$ consists of operators $M_{x}=\outerproduct*{x}{x}$, the projectors onto the computational basis states.

In the circuit model of quantum computation, we are given an input $x\in\{0,1\}^{n}$, which is a classical bit string. The first step is to prepare an $m$-qubit qubit quantum input state $|\psi\rangle$, where $m=\,$poly$(n)$. A unitary operator $U$ is then applied to $\ket*{\psi}$, and finally the output state is (customarily) measured in the computational basis – without loss of generality. The measurement outcome corresponds to a classical bit string $y\in\{0,1\}^{m}$, which is obtained with probability $\absolutevalue{\bra*{y}U\ket*{\psi}}^{2}$ and which we refer to as the output of the computation.

In practice, a quantum computer will be built using a finite set of quantum gates which act on a finite number of qubits. Typically, we consider gates acting on either one or two qubits, leaving the others invariant. A set of quantum gates is said to be universal if any unitary operator can be approximated ‘well-enough’ using only gates from this set. More precisely, a set of gates $S$ is universal if any unitary operator $U$ can be decomposed into the sequence $U_{L}U_{L-1}\ldots U_{1}$, such that $||U-U_{L}\ldots U_{1}||_{2}\leq\varepsilon$, for any $\varepsilon>0$, where the $U_{k}\in S$. There are many such universal gate sets, such as for instance the Toffoli gate (which acts on three bit/qubits) and the Hadamard gate, or single-qubit rotations with a CNOT. Thus, any arbitrary unitary operator $U$ can be implemented given a universal set of gates.

This thus tells us that any arbitrary unitary operator $U$ can be approximated by the sequence $U_{L}U_{L-1}\ldots U_{1}$ to accuracy $\varepsilon$. But, how many gates $L$ are required to achieve a good accuracy? The Solovay-Kitaev theorem (see [NC02, Appendix 3]) states that L=\mathcal{O}\big{(}\log^{2}\frac{1}{\varepsilon}\big{)}, and thus exponential accuracy can be achieved using only a polynomial number of gates.

Finally, we discuss an important tool used in quantum computation, the oracle. Here, we are given a boolean function $f:\{0,1\}^{n}\to\{0,1\}^{m}$. The function is said to be queried via an oracle $\mathcal{O}_{f}$, if given the input $|x\rangle|q\rangle$ (where $x\in\{0,1\}^{n}$ and $q\in\{0,1\}^{m}$), we can prepare the output $|x\rangle|q\oplus f(x)\rangle$, where $\oplus$ denotes addition modulo $2$. That is, the mapping

[TABLE]

can be implemented by a unitary circuit $U_{f}$, which takes the form

[TABLE]

The effect of the oracle needs to be determined on all basis states, and the definition will always be given in terms of a state $|q\rangle$.

1.4 Quantum algorithms and machine learning

Quantum algorithms, in some cases, have the capacity to achieve significant speed-ups compared to classical algorithms. Most notably, classically, the prime factorization of an $n$-bit integer using the general number field sieve is $\Omega(N^{1/3}\log^{2/3}N)$, where $N=\mathcal{O}(2^{n})$, and thus takes time exponential in the number of bits. In contrast, Shor’s factorization algorithm [Sho99] achieves an astonishing exponential speed-up with a polynomial runtime of $\mathcal{O}((\log n)^{3})$. Another impressive result is the quadratic speed-up obtained by Grover’s algorithm for unstructured database search [Gro96, Gro97], which we discuss in detail in section 2.7. These are just some examples of the many quantum algorithms which have been devised over the past decades [Cle+98, Mon15].

Machine learning [AMMIL12] has had and continues to have a significant impact for artificial intelligence and more generally for the advancement of technology. Naturally, this raises the question of whether quantum computing could enhance and improve current results and techniques. The HHL algorithm [HHL09] considers the problem of solving a system of linear equations, which on a classical computer takes time polynomial in the system size $n$. At its heart, the problem reduces to matrix inversion, and the HHL algorithm offers an exponential speed-up for this task, with a certain number of important caveats. This in turn raised the question of whether quantum algorithms could accelerate machine learning tasks, which is referred to as quantum machine learning (QML) – see the following reviews [ABG06, Cil+17].

An interesting example which illustrates how quantum computing might help with machine learning tasks is quantum recommendation systems [KP16]. Here, the netflix problem [KBV09, BK07] is considered, whereby we are given $m$ users and $n$ films, and the goal is to recommend a film to a user which they have not watched, that they would rate highly given their previous rating history. The users’ preferences can be represented by a matrix $P$ of dimension $m\times n$, where the $(i,j)^{\text{th}}$ entry corresponds to the $i^{\text{th}}$ user’s rating of the $j^{\text{th}}$ film. Of course, the elements of $P$ are not all known, and the question is to provide good recommendations to the user. Classical algorithms for this problem run in time polynomial with matrix dimension, that is \mathcal{O}\big{(}\text{poly}(MN)\big{)}. In contrast, there exists a quantum algorithm with runtime complexity scaling as \mathcal{O}\big{(}\text{poly}\log(MN)\big{)}, thus providing an exponential speed-up [KP16].

Another important example is the classical perceptron model, where we are given $N$ labeled data points which are linearly separable and the goal is to find the separating hyperplane. Classically, we have that the number of training updates scales as \mathcal{O}\big{(}\frac{1}{\gamma^{2}}\big{)}, where $\gamma$ is the margin, i.e. the shortest distance between a training point and the separating hyperplane. In contrast, the quantum perceptron model [WKS16] exploits Grover’s search algorithm (see section 2.7), in order to achieve a quadratic speed-up \mathcal{O}\big{(}\frac{1}{\gamma}\big{)}.

Another important classical model for supervised learning is support vector machines (SVM), which are used for data classification and regression. For a special case of SVMs, known as least-squares SVMs, the classical runtime is polynomial in the number of training examples $N$ and their dimension $d$, $\mathcal{O}(\log\frac{1}{\varepsilon}\text{poly}(d,N))$, where $\varepsilon$ denotes the accuracy. In contrast, quantum SVM [RML14] offer an exponential speed-up in the dimensions and input number with a runtime of $\mathcal{O}(\frac{1}{\epsilon}\log(dN))$.

Finally, we note that quantum algorithms have also been developed for unsupervised learning, that is, in the case of unlabeled data [ABG13, WKS14], which also present a speed-up over their classical counterparts.

All of the algorithms mentioned here use the HHL – or some other quantum linear systems algorithm – as a subroutine.

1.5 Structure

We have aimed for a mostly modular presentation, so that a reader familiar with a particular subroutine of the HHL algorithm, say, phase estimation, can skip the relevant section if they wish. The structure of the text goes as follows.

First, in section 2, we review some of the key components used in quantum algorithms, namely notation conventions 2.1, the quantum Fourier transform 2.2, Hamiltonian simulation 2.3, quantum phase estimation 2.5, phase kickback 2.6, amplitude amplification 2.7, the uncompute trick 2.8 and finally quantum RAM 2.9. Next, in section 3, we present a detailed discussion of the HHL algorithm. We first formally define the problem in 3.1, before then discussing the algorithm in detail in 3.2, along with an error analysis in 3.3. In section 3.5 we consider the computational complexity of the problem, and in 3.5 its optimality. Then, in 3.6, the algorithm is extended to the case of non-Hermitian matrices.

Next, in section 4, we introduce modern updates to the algorithm, namely: improvements in the dependency on the condition number in 4.1; improvements on the precision number in 4.2 and in section 4.3, further improvements based on quantum singular value estimation giving an algorithm for dense matrices. More specifically, we discuss Jordan’s lemma and its consequences in LABEL:subsec:intersection_subspaces, before reviewing the singular value decomposition in 4.3.1 and finally presenting the quantum singular value estimation 4.3.2 and its application to linear systems in 4.3.3. This section deviates from the pedagogical style of the previous sections, giving an overview of the important ideas behind the results as opposed to all of the gory details.

2 Quantum algorithms: fundamental components

Here we review some of the fundamental ideas routinely used in the design of quantum algorithms, which will subsequently be applied in the HHL algorithm.

2.1 Notation and conventions

Any integer $k$ between [math] and $N-1$, where $N=2^{n}$ may be expressed as an $n$-bit string $k=k_{1}\ldots k_{n}$, i.e. $k=\sum_{l=1}^{n}k_{l}2^{n-l}$. Furthermore, given an integer $j$, it is easy to verify that $\frac{j}{2^{m}}=0.j_{n-m+1}\ldots j_{n}$.

The Hadamard matrix is defined as H:=\frac{1}{\sqrt{2}}\smallmatrixquantity[1&1\\ 1&-1]. Given a vector $\textbf{x}=(x_{1},\ldots,x_{n})$, the Euclidean norm of x is $\norm{\textbf{x}}_{2}=\sum\sqrt{|z_{1}|^{2}+\ldots+|z_{n}|^{2}}$. For an $m\times n$ matrix $A$ with elements $a_{ij}$, the operator norm is given by $\norm{A}_{2}=\sum_{j=1}^{n}\sqrt{\sum_{i=1}^{m}|a_{ij}|}$, i.e., it is the sum of the Euclidean norms of the column vectors of $A$. Note that in both cases this is equivalent to $\ell_{2}$-norm. Next, we present the QFT.

2.2 Quantum Fourier transform

The QFT is at the heart of many quantum algorithms. It is the quantum analogue of the discrete Fourier transform, see 2.2.1, and is presented in section 2.2.2. In section 2.2.3, we will see how the QFT may be implemented efficiently using a quantum computer.

2.2.1 The discrete Fourier transform

The Fourier transform is an important tool in classical physics and computer science, as it allows for a signal to be decomposed into its fundamental components, i.e. frequencies. The Fourier transform tells us what frequencies are present and to what degree.

In the discrete setting, we have that the DFT is a square invertible matrix $D$ of dimension $N$, where $D_{jk}=\frac{1}{\sqrt{N}}\omega^{jk}$ and $\omega=e^{i2\pi/N}$. It is easy to show that the columns of this matrix are orthogonal and have unit length, and thus the set of column vectors form an orthonormal basis which we refer to as the Fourier basis. If the DFT is applied to a vector using matrix multiplication, then the time complexity scales as $\mathcal{O}(N^{2})$. Crucially, by using a divide and conquer approach this can be improved to $\mathcal{O}(N\log N)$, which is referred to as the fast Fourier transform (FFT).

2.2.2 The quantum Fourier transform

In a similar way, the QFT is defined by mapping each computational basis state $|j\rangle$ to a new quantum state |f_{j}\rangle=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}\exp\Big(i\frac{2\pi jk}{N}\Big{missing})|k\rangle. The set of orthonormal states $\{|f_{j}\rangle\}$ form an orthonormal basis set called the Fourier basis. The quantum Fourier transform with respect to an orthonormal basis $\ket*{x}\in\{\ket*{0},\ldots,\ket*{N-1}\}$ is defined as the linear operator with the following action on the basis vectors:

[TABLE]

The inverse Fourier transform is then defined as

[TABLE]

But, what does this mean in terms of individual qubits? First, we represent the integer $k$ in binary notation, $k=\sum_{l=1}^{n}k_{l}2^{n-l}$, and thus the Fourier basis states can be expressed as: |f_{j}\rangle=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}\exp\Big(i2\pi j\sum_{l=1}^{n}k_{l}2^{-l}\Big{missing})|k\rangle. Expanding this expression, we have

[TABLE]

Expanding the summation gives

[TABLE]

Finally, the operation $\frac{j}{2^{m}}=0.j_{n-m+1}\ldots j_{n}$ corresponds to the decimal expansion of $j$ up to $m$ bits, and we can thus write

[TABLE]

We initially applied the QFT to the computational basis state $|j\rangle=|j_{1}\ldots j_{n}\rangle$. From Eq. (10), we see that information pertaining to the input $j$ is disseminated throughout the relative phase on each individual qubit. Thus, given the final state $|f_{j}\rangle$, applying the inverse QFT would yield the input string $j$. Equivalently, one could obtain $j$ by performing a measurement in the Fourier basis.

2.2.3 Implementation of the QFT

The goal is to obtain the quantum state $|f_{j}\rangle$ after applying a quantum circuit to an all-zero input state. From Eq. (10) we see that the state $|f_{j}\rangle$ corresponds to a state where each qubit is initialised in the state $|+\rangle=\frac{|0\rangle+|1\rangle}{\sqrt{2}}$ and subsequently acquires a relative phase of \exp\big(i2\pi j2^{-k}\big{missing}), where $1\leq k\leq n$, and where we recall that $1\leq j\leq N-1$.

We now consider the quantum circuit which can implement this state. First, it is easy to see that a state of the form \frac{1}{\sqrt{2}}\Big{(}|0\rangle+\exp\big(i2\pi 0.j_{n}\big{missing})|1\rangle\Big{)} corresponds to either the state $|+\rangle$ or $|-\rangle$ depending on the value of $j_{n}\in\{0,1\}$. This can be expressed as $H|j_{n}\rangle$, and can thus be obtained by the application of a Hadamard gate $H$ to a qubit in the state $|j_{n}\rangle$. Next, the state of the second qubit is given by \frac{1}{\sqrt{2}}\Big{(}|0\rangle+\exp\big(i2\pi 0.j_{n-1}j_{n}\big{missing})|1\rangle\Big{)}, which can be re-expressed as \frac{1}{\sqrt{2}}\Big{(}|0\rangle+(-1)^{j_{n-1}}\exp\big(i2\pi\frac{j_{n}}{2^{2}}\big{missing})|1\rangle\Big{)}. Thus, this corresponds to first preparing the state $H|j_{n-1}\rangle$, and then applying a controlled rotation to the qubit, where the control is the $n$th qubit in the state $|j_{n}\rangle$. Thus, the state on the first two qubits can be obtained by preparing the state $|j_{n-1}j_{n}\rangle$, applying the Hadamard gate to the $n$th qubit, and then a controlled rotation $R_{2}$ with qubit $|j_{n-1}\rangle$ as control, where we have that $R_{k}$ is given by:

[TABLE]

and controlled by qubit $|j_{n-k+1}\rangle$. Finally, SWAP operations are performed throughout for the qubits to be in the correct order. This approach can be extended to the $n$ qubits, where the number of gates scales as $\mathcal{O}(n^{2})$, and thus the QFT can be efficiently implemented in a quantum circuit.

2.3 Hamiltonian simulation

Most quantum algorithms for machine learning, and in particular the HHL algorithm, leverage quantum Hamiltonian simulation as a subroutine. Here, we are given a Hamiltonian operator $\hat{H}$, which is a Hermitian matrix, and the goal is to determine a quantum circuit which implements the unitary operator $U=e^{-i\hat{H}t}$, up to given error. The evolution of a quantum state $\ket*{\Psi}$ under a unitary operator is given, for simplicity, by the time-independent Schrödinger equation:

[TABLE]

the solution to which can be written as $\ket*{\Psi(t)}=\exp(-i\hat{H}t)\ket*{\Psi}$.

Depending on the input state and resources at hand, there exists many different techniques to achieve this [BC09, BCK15, LC17, LC17a, LC16]. We give a brief introduction to this large field of still ongoing research, and interested readers can find further details in the lecture notes of Childs [Chi17][Chapter V] and the seminal work [Chi+03].

The challenge is due to the fact that the application of matrix exponentials are computationally expensive. For instance, naive methods require time $\mathcal{O}\left(N^{3}\right)$ for a $N\times N$ matrix, which is restrictive even in the case of small size matrices. In the quantum case, the dimension of the Hilbert space grows exponentially with the number of qubits, and thus any operator will be of exponential dimension. Applying such expensive matrix exponentials has been studied classically, and in particular the classical simulation of such time-evolutions is known to be hard for generic Hamiltonians $\hat{H}$. As a consequence, new more efficient methods need to be introduced. In particular, a quantum computer can be used to simulate the Hamiltonian operator, a task known as Hamiltonian simulation, which we wish to perform efficiently. More specifically, we can now define an efficient quantum simulation as follows:

Definition 1**.**

(Hamiltonian Simulation)* We say that a Hamiltonian $\hat{H}$ that acts on $n$ qubits can be efficiently simulated if for any $t>0,\epsilon>0$, there exists a quantum circuit $U_{\hat{H}}$ consisting of $\text{poly}(n,t,1/\epsilon)$ gates such that $\norm{U_{\hat{H}}-e^{-i\hat{H}t}}<\epsilon$. Since any quantum computation can be implemented by a sequence of Hamiltonian simulations, simulating Hamiltonians in general is $\mathsf{BQP}$-hard, where $\mathsf{BQP}$ refers to the complexity class of decision problems efficiently solvable on a universal quantum computer [KSV02].*

Note that the dependency on $t$ is important and it can be shown that at least time $\Omega(t)$ is required to simulate $\hat{H}$ for time $t$, which is stated formally by the no fast-forwarding theorem [Ber+07]. There are, however, no nontrivial lower bounds on the error dependency $\epsilon$. The hope to simulate an arbitrary Hamiltonian efficiently is diminished, since it NP-hard to find an approximate decomposition into elementary single- and two-qubit gates for a generic unitary and hence also for the evolution operator [SBM06, Ite+16, Kni95]. Even more so, the optimal circuit sythesis was even shown to be QMA-complete [JWB03]. However, we can still simulate efficiently certain classes of Hamiltonians, i.e. Hamiltonians with a particular structure. One such example is the case when $\hat{H}$ only acts nontrivially on a constant number of qubits, as any unitary evolution on a constant number of qubits can be approximated with error at most $\epsilon$ using $\text{poly}(\log(1/\epsilon)$ one- and two-qubit gates, on the basis of Solovay-Kitaev’s theorem. The origin of the hardness of Hamiltonian simulation stems from the fact that we need to find a decomposition of the unitary operator in terms of elementary gates, which in turn can be very hard for generic Hamiltonians. If $\hat{H}$ can be efficiently simulated, then so can $c\hat{H}$ for any $c=\text{poly}(n)$ [Chi17]. In addition, since any computation is reversible, $e^{-i\hat{H}t}$ is also efficiently simulatable and this must hence also hold for $c<0$.

Finally, we note that the definition of efficiently simulatable Hamiltonians further extends to unitary matrices, since every operator $U_{\hat{H}}$ corresponds to a unitary operator, and furthermore every unitary operator can be written in the form $\exp(i\hat{H})$ for a Hermitian matrix $\hat{H}$. Hence, we can similarly speak of an efficiently simulatable unitary operator, which we will use in the following.

2.3.1 Trotter-Suzuki methods

For any efficiently simulatable unitary operator $U$, we can always simulate the Hamiltonian $\hat{H}$ in a transformed basis $U\hat{H}U^{\dagger}$, since

[TABLE]

which follows from the fact that if $U=U^{\dagger}$ is unitary, then we have that $(U\hat{H}U^{\dagger})^{m}=U\hat{H}^{m}U^{\dagger}$, which can easily be proven by induction. Another simple but useful trick is given by the fact that, given efficient access to any diagonal element of a Hamiltonian $H_{ii}=\bra*{i}\hat{H}\ket*{i}$, we can simulate the diagonal Hamiltonian using the following sequence of operations. Let $\Rightarrow$ indicate a computational step, such that we can denote in the following a sequence of maps to a state:

[TABLE]

In words, we first load the entry $H_{ii}$ into the second register, then apply a conditional-phase gate $\exp(-i{H}_{ii}t)$ and then reverse the loading procedure to set the last qubit to zero again. Since we can apply this to a superposition and using linearity, we can simulate any efficiently diagonalisable Hamiltonian. More generally, any $k$-local Hamiltonian, i.e. a sum of polynomially many terms in the number of qubits that each act on at most $k=\mathcal{O}\left(1\right)$ qubits, can be simulated efficiently. Indeed, since each of the terms in the sun acts only on a constant number of qubits, it can be efficiently diagonalised and thus simulated. In general, for any two Hamiltonian operators $\hat{H}_{1}$ and $\hat{H}_{2}$ that can be efficiently simulated, the sum of both $\hat{H}_{1}+\hat{H}_{2}$ can also be efficiently simulated, as we will argue below, first for the commuting case and then for the non-commuting case.

This is trivial if the two Hamiltonians commute. Indeed, we now omit the coefficient $it$ and consider for simplicity the operator $\exp(\hat{H}_{1}+\hat{H}_{2})$. By applying a Taylor expansion, followed by the Binomial theorem and the Cauchy product formula (for the product of two infinite series), we have

[TABLE]

Note that this is only possible since for the Cauchy formula we can arrange the two terms accordingly and do not obtain commutator terms in it. However (recall the famous Baker-Campbell-Hausdorff formula, see e.g. [Ros02]), this is not so for the general case, i.e. when the operators don’t commute. Here we need to use the Lie-Product formula [Ros02]:

[TABLE]

If we want to restrict the simulation to a certain error $\epsilon$, it is sufficient to truncate the above product formula after a certain number of iterations $m$, which we will call the number of steps:

[TABLE]

which, as we will show, can be achieved by taking $m=\mathcal{O}\left((\max{(||\hat{H}_{1}||,||\hat{H}_{2}||)})^{2}t^{2}/\epsilon\right)$, where we require that $\max{(||\hat{H}_{1}||,||\hat{H}_{2}||)}\sim\mathcal{O}\left(\text{poly}(n)\right)$ for the evolution to be efficiently simulable. To see this, observe that, from the Taylor expansion,

[TABLE]

We need to expand a product of the form $(A+B)^{m}$, where the operators $A$ and $B$ are non-commuting. Thus, we have that

[TABLE]

where $A$ and $B$ do not commute in general. Specifically we have $A=e^{-i(\hat{H}_{1}+\hat{H}_{2})t/m}$ and $B=\mathcal{O}\left(\frac{t^{2}\max{\left(\hat{H}_{1},\hat{H}_{2}\right)}^{2}}{m^{2}}\right),$ and the $m$ first order terms in $B$ have the form

[TABLE]

for $k\in[0,m-1]$. Next, we consider terms of order greater than one, where we have powers of $B^{l}$ for $l>2$. Let us note that for $m>1$, we have $C/m^{2}+C/m^{3}+C/m^{4}+\ldots\leq\mathcal{O}\left(C/m^{2}\right)$ for $m>1$. Thus, these first order and greater terms can be absorbed in the $\tilde{\mathcal{O}}$ notation. Furthermore, in the following, we do not explicitly write the exponentials of the form $\exp(-i(H_{1}+H_{2})t(m-k-1)/m)$ as these will not play a role for bounding the error in the norm due to their unitarity. So, we continue to bound

[TABLE]

We can now finally consider the error of the simulation scheme, see Eq. (22), which using Eq. (27), yields

[TABLE]

In order to have this error less than $\epsilon$, the number of steps $m$ must be $m=\mathcal{O}\left(\frac{t^{2}\max{(\hat{H}_{1},\hat{H}_{2})}^{2}}{\epsilon}\right)$.

This is a naive and non-optimal scheme. It can be shown that one can use higher-order approximation schemes, such that $\hat{H}_{1}+\hat{H}_{2}$ can be simulated for time $t$ in $t^{1+\delta}$ for any positive but arbitrarily small $\delta$ [Ber+07, BC09].

These so-called Trotter-Suzuki schemes can be generalized to an arbitrary sum of Hamiltonians which then leads to an approximation formula given by

[TABLE]

The following definitions are useful:

Definition 2**.**

(Sparsity)* An $N\times N$ matrix is said to be $s$-sparse if it has at most $s$ entries per row.*

Definition 3**.**

(Sparse matrix)* An $N\times N$ matrix is said to be sparse if it has at most $\operatorname{poly}\left(\log N\right)$ entries per row.*

Note that the sparsity depends on the basis in which the matrix is given. However, given an arbitrary matrix, we do not a priori know the basis which diagonalises it (and hence gives us a diagonal matrix with sparsity $1$), and so we need to deal with a potentially dense matrix, i.e. a matrix which has $N$ entries per row.

Definition 4**.**

(Row computability)* The entries of a matrix $A$ are efficiently row computable if, given the indices $i,j$, we can obtain the entries $A_{ij}$ efficiently, i.e. in $\mathcal{O}\left(s\right)$ time, where $s$ is the sparsity as defined above.*

2.3.2 Graph-colouring method

Crucially, the simulation techniques described above can allow us to efficiently simulate sparse Hamiltonians. Indeed, if for any index $i$, we can efficiently determine all of the indices $j$ for which the term $\bra*{i}\hat{H}\ket*{j}$ is nonzero, and furthermore efficiently obtain the values of the corresponding matrix elements, then we can simulate the Hamiltonian $\hat{H}$ efficiently, as we will describe below.

This method of Hamiltonian simulation is based on ideas from graph theory, and we will now first briefly introduce a couple of key notions relevant in our discussion. For further information, we refer the reader to the existing literature [Wes+01]. An undirected graph $G=(V,E)$ is specified by a set $V$ of $|V|=N$ vertices and a set $E$ of $|E|=M$ edges, i.e., unordered pairs of vertices. When two vertices form an edge they are said to be connected. A graph can be represented by its adjacency matrix $A$, where $A_{ij}=1$ if the vertices $i$ and $j$ are connected, and $A_{ij}=0$, otherwise. The degree of a vertex is given by the number of vertices it is connected to. The maximum degree of a graph refers the maximum degree taken over the set of vertices. The problem of edge colouring considers if, given $k$ colours, each edge can be assigned a specific colour with the requirement that no two edges sharing a common vertex should be assigned the same colour. Vizing’s theorem tells us that, for a graph with maximum degree $d$, an edge colouring exists with $k$ at most $d+1$. Finally, a bipartite graph is a graph, where the set of vertices can be separated into two disjoint subsets $V_{1}$ and $V_{2}$ such that $V=V_{1}\cup V_{2}$ and no two vertices belonging to the same subset are connected, i.e., for every $i,j\in V_{k}$, we have $(i,j)\notin V_{k}$, for $k=1,2$.

Previously, we saw that a Hamiltonian operator can be represented by a square matrix. Thus, a graph can be associated with any Hamiltonian by considering the adjacency matrix with a $1$ at every non-zero entry of the Hamiltonian, and a [math] elsewhere, in the spirit of combinatorial matrix theory. For a matrix of dimension $N$, this will thus correspond to a graph with $N$ vertices. Previously, we saw that sparse Hamiltonians have at most polylog$(N)$ entries per row, and thus $\mathcal{O}\left(n\right)$, for $n=\log N$ entries in total. This will translate into a graph having a number of edges $|E|=\mathcal{O}\left(n\right)$.

Childs [Chi+03] proposed an efficient implementation for the simulation of sparse Hamiltonians by using the Trotterization scheme presented above (c.f. section 2.3.1) and a local colouring algorithm of the graph associated with the $s$-sparse Hamiltonian. The core idea is to colour the edges of the Hamiltonian $\hat{H}$. Then, the Hamiltonians corresponding to each subgraphs defined by a specific colour can be simulated, and finally the original Hamiltonian recovered via the Trotter-Suzuki method [Chi+03].

More precisely, the first step is to find an edge-colouring of the graph associated with the Hamiltonian $\hat{H}$. This will be achieved using $k$ colours, which in the case of a sparse Hamiltonian, will be at most polylog$|n|$. Next, the graph can be decomposed by considering the subgraphs corresponding to a single colour. We thus obtain a decomposition of the original Hamiltonian in a sum of sparse Hamiltonians, containing at most polylog$|n|$ terms. It is easy to convince oneself that each of these terms consists of a direct sum of two-dimensional blocks. Indeed, each adjacency matrix corresponding to a subgraph will be symmetric, with at most one entry per row, meaning that the evolution on any one of these subgraphs takes place in isolated two-dimensional subspaces. Thus, each Hamiltonian term can be easily diagonalised and simulated using the diagonal Hamiltonian simulation procedure as given in Eq. (13).

A crucial step in this procedure is the classical algorithm for determining the edge colouring efficiently. Vizing’s theorem guarantees the existence of an edge colouring using $d+1$ colours. But, the question remains as to how this can be efficiently achieved. Indeed, even though we are given the adjacency matrix representation of the entire graph, we will now restrict ourselves to accessing only local information i.e. each vertex has only access to information regarding it nearest-neighbours. Finding an optimal colouring is an NP-complete problem. However, there are polynomial time algorithms that construct optimal colourings of bipartite graphs, and colourings of non-bipartite simple graphs that use at most $d+1$ colours. It is important to note that the general problem of finding an optimal edge colouring is NP-hard and the fastest known algorithms for it take exponential time.

We thus now present a local edge-colouring scheme achieving a $d^{2}$-colouring (where we recall that $d$ is the maximum degree of the graph, i.e. the sparsity) for the case of a bipartite graph. This, using a reduction [Chi17], is sufficient for the simulation of an arbitrary Hamiltonian. Crucially, this scheme is efficient if the graph is sparse, i.e. if $d=$polylog$(n)$. We note that better schemes exist [Ber+07, BC09, BCK15] and can allow for polynomial improvements of the simulation scheme in comparison to the one given here.

Lemma 1**.**

(Efficient bipartite graph colouring [Lin87, Lin92])*

Suppose we are given an undirected, bipartite graph $G$ with $n$ vertices and maximum degree $d$ (i.e. each vertex is connected to a maximum of $d$ other vertices - the so called neighbours - which is similar to sparsity $s$), and that we can efficiently compute the neighbours of any given vertex. Then there is an efficiently computable edge colouring of $G$ with at most $d^{2}$ colours.*

Proof.

The vertices of $G$ are ordered and numbered from $1$ through $n$. For any vertex $a$, let $index(a,b)$ denote the index of vertex $b$ in the list of neighbours of $a$, ordered in increasing number. For example, let $a$ have the neighbours $c,d$ with the list of neighbours neighbours $(a):=\{c,d\}$ of $a$. Then, we have that index$(a,c)=1$, and index$(a,d)=2$. Then define the colour of the edge $(ab)$, where $a$ is from the left part of the bipartition and $b$ is from the right for all $a$ and $b$ which have an edge. The colouring of this edge $(ab)$ is then assigned to be the ordered pair colour$(ab):=($index$(a,b),$ index$(b,a))$. Recall that an edge colouring assigns a colour to each edge so that no two adjacent edges share the same colour. These assigned colours in form of the index-pairs give a valid colouring since if $(a,b)$ and $(a,d)$ have the same colour, then index$(a,b)=\,$index$(a,d)$, so $b=d$ and similarly, if $(a,b)$ and $(c,b)$ have the same colour, then index$(b,a)=$ index$(b,c)$, so $a=c$. ∎

Using this lemma we can then perform Hamiltonian simulation in the following manner. First we ensure that the associated graph is bipartite by simulating the evolution according to the Hamiltonian $\hat{B}=\sigma_{x}\otimes\hat{H}$, a block-anti-diagonal matrix

[TABLE]

The graph associated to this will be bipartite, with the same sparsity as $\hat{H}$ [Chi17]. Observe that simulating this reduces to simulating $\hat{H}$ since

[TABLE]

Without loss of generality let us assume now that $\hat{H}$ has only off diagonal entries. Indeed, any Hamiltonian $\hat{H}$ can be decomposed as the sum of diagonal and off-diagonal terms $\hat{H}_{\mathrm{diagonal}}+\hat{H}_{\mathrm{off-diagonal}}$, which can then be simulated the sum using the rule given in Eq. (21). We can then, for a specific vertex $x$ and a colour $c$, compute the evolution by applying the following three steps:

First we compute the complete list of neighbours of $x$ (i.e. the neighbour list and the indices) and each of the colours associated to the edges connecting $x$ with its neighbours, using the above local algorithm for graph colouring from Lemma 1. 2. 2.

Let $v_{c}(x)$ denote the vertex adjacent to $x$ via an edge with colour $c$. We then, for a given $x$, compute $v_{c}(x)$ and retrieve the Hamiltonian matrix entry $H_{x,v_{c}(x)}$. We can then implement the following quantum state $\ket*{x,v_{c}(x),H_{x,v_{c}(x)}}$ i.e. three qubit registers in which we load the elements $x$ in the first one, $v_{c}(x)$ in the second and then load the matrix element into the last one. More specifically, we here prepare the quantum state $\ket*{x,v_{c}(x)}\otimes\ket*{a}$ and then using the state preparation oracle (e.g. qRAM, see section 2.9), we obtain $\ket*{x,v_{c}(x)}\otimes\ket*{a\oplus H_{x,v_{c}(x)}}$ which can be done efficiently, i.e. in time $\mathcal{O}\left(\log(n,m)\right)$, which is the time required to access the data. 3. 3.

We then simulate the ($\hat{H}$-independent, i.e. only depending on the local entry of $\hat{H}$ but not the general matrix) Hamiltonians $\hat{h}$, i.e. the at most polylog$|n|$ Hamiltonians we obtain from the graph colouring. Note that this is a sparse Hamiltonian which acts only on the $d$ neighbours as described in the colouring step. The simulation is efficient, since $\hat{h}$ can be diagonalised in constant time, as it consists of a direct sum of two-dimensional blocks. Next, we apply a scheme, described below, whereby each complex matrix entry is decomposed into a real part $x$ and imaginary part $y$ and simulated separately. We can then simulate the diagonalised Hamiltonians such that we implement the following mapping

[TABLE]

where $h_{a}$ is a diagonal element of the diagonalised Hamiltonian $\hat{h}$. This can also be done in superposition.

The Hamiltonian to be simulated has complex entries, and can thus be decomposed in real and imaginary parts. Let $v_{c}(a)$ denote the vertex connected to $a$ via an edge of colour $c$. The original Hamiltonian had complex entries, and we can express the entry associated to vertex $v_{c}(a)$ as a sum of a real part $x_{c}(a)$ and imaginary part $y_{c}(a)$. If we assume that these can be loaded independently i.e. $H_{a,v_{c}(a)}=x_{c}(a)+i\ y_{c}(a)$ can be loaded separately, then we can introduce the oracles $V_{c},W_{c}$ which allow for the following mappings to be implemented:

[TABLE]

and similarly the inverse operations, for which it holds that $V_{c}^{\dagger}=V_{c}$, $W_{c}^{\dagger}=W_{c}$, since we have bitwise adding modulo $2$ in the loading procedure.

In order to simulate the complex entries we need to implement a procedure which allows us to apply both parts individually and still end up in the same basis-element such that these sum up to the actual complex entry, i.e. that we can apply $x_{c}(a)$ and $i\ y_{c}(a)$ separately to the same basis state. We use multiple steps to do so.

Given the above oracles, we can similarly simulate the following Hermitian operations (note that this is not a unitary operation):

[TABLE]

where we apply the Hermitian operators to multiply with $x$ (and $iy$) and the swap operation to the registers, which can be implemented efficiently since the swap can be done efficiently. The operator $T$ is described in detail in [Chi+03] We then can implement the operator

[TABLE]

where the sum is about all colours $c$. This acts then on $\ket*{a,0,0}$ as $H$, since

[TABLE]

which can be confirmed using the fact that $v_{c}(v_{c}(a))=a$ and since we have modulo $2$ addition, i.e. $a\oplus a=0$.

2.4 Erroneous Hamiltonian simulation

One might wonder what would happen with errors in the Hamiltonian simulator. For example, imagine that we simulate the target Hamiltonian with simulator such as a quantum computer, and this simulator introduces some random error terms. This will be an issue for as long as we do not have fully error corrected quantum computers, or if we use methods like quantum density matrix exponentiation which can have errors in the preparation of the state we want to exponentiate (see [LMR13, Kim+17]). For example, this is relevant for a method called sample-based Hamiltonian simulation [LMR13, Kim+17] where we perform the quantum simulation of a density matrix $\rho$, i.e. trace-$1$ Hermitian matrix, which can have some errors.

For a more in depth introduction and analysis see for example [CMP17]. Errors in the Hamiltonian simulator have been investigated in depth and here, we only want to give the reader some tools to grasp how one could approach such a problem in this small section. For more elaborate work on this we refer the reader to [CMP17].

In the following we will use so called matrix Chernoff-type bounds. Let us recall some results from statistics:

Theorem 1** (Bernstein [Ber27]).**

Let $\{X_{i}\}_{i=1}^{n}$ be random variables and let $X=\sum_{i}X_{i}$, such that $\mathbb{E}(X)$ is the expectation value and $\mathbb{E}(X^{2})=\sigma^{2}$, $|X|\leq M$ is bounded and $X_{i}$ are independent, identically distributed copies. Then, for all $t>0$,

[TABLE]

This fundamental theorem in statistics is making use of the independence of the sampling process in order to obtain a concentration of the result in high probability. In order to provide bounds for Hamiltonian simulation, we will need to use matrix versions of these Chernoff-style results. We will thereby make certain assumptions about the matrix, such as for instance that it is bounded in norm, and that the matrix variance statistic - a quantity that is a generalization of the variance - has a certain value. We now state first the result and then prove it.

Lemma 2**.**

(Faulty Hamiltonian simulator)* Hamiltonian simulation of a $N\times N$ Hamiltonian operator $\hat{H}=\sum_{i}\hat{H}_{i}$ with a faulty simulator that induces random error terms $\{\hat{H}_{i}^{err}\}$ (random matrices) with expectation value $\mathbb{E}[\hat{H}^{err}_{i}]=0$ with bounded norm $\norm{\hat{H}^{err}_{i}}_{2}\leq R$ for all $i$ and bounded matrix variance statistic $v(\hat{H}^{err}:=\sum_{i}\hat{H}^{err}_{i})=\max{\left(\norm{\mathbb{E}(\hat{H}^{err}\hat{H}^{err\dagger})}_{2},\norm{\mathbb{E}(\hat{H}^{err\dagger}\hat{H}^{err})}_{2}\right)}$ in each term of the simulation can be simulated with an error less than $\mathcal{O}((1+t)\cdot\epsilon)$ using $m=\mathcal{O}\left(\max_{i}\norm{\hat{H}_{i}}_{2}^{2}t^{2}/\epsilon\right)$ steps with probability of at least*

[TABLE]

Proof.

To prove this we will need a theorem that was developed independently in the two papers [Tro12, Oli09], which is a matrix extension of Bernstein’s inequality. Recall that the standard Bernstein inequality is a concentration bound for random variables with bounded norm, i.e. it tells us that a sum of random variables remains, with high probability, close to the expectation value. This can be extended to matrices which are drawn form a certain distribution and have a given upper bound to the norm. We call a random matrix independent and centered if each entry is independently drawn from the other entries and the expectation of the matrix is the zero matrix.

Theorem 2** (Matrix Bernstein).**

Let $X_{1},...,X_{n}$ be independent, centered random matrices with dimensionality $N\times N$, and assume that each one is uniformly bounded

[TABLE]

We define the sum of these random matrices $X=\sum_{i=1}^{n}X_{i}$, with matrix variance statistic of the sum being

[TABLE]

where $\mathbb{E}(XX^{\dagger})=\sum_{i=1}^{n}\mathbb{E}(X_{i}X_{i}^{\dagger})$. Then

[TABLE]

Furthermore,

[TABLE]

Let us then recall the error in the Hamiltonian simulation scheme from above.

[TABLE]

We then use the assumption that the matrix variance statistic is bounded and that all the $H^{err}_{i}$ are bounded in norm by $R$ in order to be able to apply the above theorem. Using theorem 2 in order to probabilistically bound the first term, and observing that for $m=\mathcal{O}\left(\max_{i}\norm{H_{i}}^{2}t^{2}/\epsilon\right)$ we can bound the second term by $\epsilon$, assuming that $\max_{i}{H^{err}_{i}}<\min_{i}{H_{i}}$ and using theorem 2 we achieve the proposed Lemma. ∎

For a more in depth analysis of Hamiltonian simulation and errors we refer the reader to [CMP17].

2.4.1 Modern methods for Hamiltonian simulation

Modern approaches like fractional-query model [Cle+09] are more complicated and we will not describe these here. However, these allow for tighter bounds and faster Hamiltonian simulation, with improved dependency on all parameters.

A conceptually different approach to Hamiltonian simulation based on Szegedy’s quantum walk [Sze04] uses the notion of a discrete-time quantum walk that is closely related to any given time-independent Hamiltonian and applies phase estimation in order to simulate the evolution. This approach has the best known performance as a function of the sparsity $s$ and evolution time $t$ but has a worse $\epsilon$-dependency, i.e. in the error. This improved method for Hamiltonian simulation scales as $\mathcal{O}\left(s\norm{H_{max}}t/\sqrt{\epsilon}\right)$ for a fixed Hamiltonian.

Other recent results based on different methods approach optimality, i.e. linear dependency in the parameters. The fact that such exponentials of the Hamiltonian can be easily performed on a quantum computer is essential to the HHL algorithm because it allows us to perform eigenvalue estimation, which we will discuss below.

2.5 Quantum phase estimation

The goal of quantum phase estimation [Kit95] is to obtain a good approximation of an operator’s eigenvalue given the associated eigenstate. Here, we consider a unitary operator $U$ acting on an $m$-qubit state, with a set of given eigenvectors $|\psi_{i}\rangle$ and associated unknown eigenvalues $\lambda_{i}$. For simplicity, we will consider a particular given eigenvector $|\psi\rangle$ with associated unknown eigenvalue $\lambda$. This eigenvalue is a complex number and we can thus write $\lambda=e^{i2\pi\phi}$, where $0\leq\phi\leq 1$ is referred to as the phase. Thus, we wish to determine a good $n$-bit approximation $\tilde{\phi}=0.\tilde{\phi}_{1}\ldots\tilde{\phi}_{n}$ of $\phi$, which will thus allow for a good $n$-bit approximation $\tilde{\lambda}$ of $\lambda$. This will be achieved by requiring $n$ ancillary qubits. The intuition is to encode this approximation within relative phases of the qubits, as we previously saw with the QFT, see section 2.2.

In first instance, we note that $U^{j}|\psi\rangle=\lambda^{j}|\psi\rangle=e^{i2\pi\phi j}|\psi\rangle$. To start with, we are given $n$ qubits prepared in the $|0\rangle^{\otimes n}$ state and an $m$-qubit quantum state intialised in the $|\psi\rangle$ state. Next, $n$ Hadamard gates are applied to the $n$ qubits in the first register, which results in the state $\frac{1}{\sqrt{N}}\sum_{x}|x\rangle|\psi\rangle$, where $N=2^{n}$.

Then, in order to obtain the $n$-bit approximation of the eigenvalue, we will apply $n$ unitary operators to the state $\psi\rangle$. More specifically, the controlled-$U^{2^{k}}$ operator is applied, where $k$ is the control qubit, with $k=0$ for the first up to $k=n-1$ for the last qubit, as illustrated in Figure 1. The sequence of controlled-$U^{2^{k}}$ operations used in phase estimation can be implemented efficiently using the technique of modular exponentiation, which is discussed in depth in [NC02][Box 5.2, Ch. 5.3.1]. This results in the state:

[TABLE]

In the following, we wish to consider $\phi$ to $n$ bits of accuracy, where the $n$-bit approximation of $\phi$ is given by $\sum_{s}\phi_{s}2^{-s}$. In all generality, the phase can be expressed as:

[TABLE]

where we have that $j=j_{n-1}\ldots j_{0}$ or alternatively $\frac{j}{2^{n}}=0.j_{1}\ldots j_{n}$. Substituting into 50 and dividing numerator and denominator by $2^{n}$ we have

[TABLE]

which it is easy to recognise is simply the QFT applied to the state $|\phi_{1}\cdots\phi_{n}\rangle$. Thus, performing a measurement in the Fourier basis gives the bit string $\phi_{1}\cdots\phi_{n}$. In the case where $\delta=0$, the measurement will yield the original state $|\phi_{1}\cdots\phi_{n}\rangle$. If $\delta\neq 0$, then the measurement will have to be repeated, in which case an upper bound on the necessary number of repetitions can be obtained.

Thus, quantum phase estimation relies on preparing $n$ ancillary qubits in equal superposition, applying $n$ controlled unitary operators to the given eigenstate, and finally performing a measurement in the Fourier basis, as illustrated in Figure 2.

We can summarize the phase estimation procedure in the following Theorem:

Theorem 3** (Phase estimation [Kit95]).**

Let unitary $U\ket*{v_{j}}=\exp(i\theta_{j})\ket*{v_{j}}$ with $\theta_{j}\in[-\pi,\pi]$ for $j\in[n]$. There is a quantum algorithm that transforms $\sum_{j\in[n]}\alpha_{j}\ket*{v_{j}}\mapsto\sum_{j\in[n]}\alpha_{j}\ket*{v_{j}}\ket*{\tilde{\theta}_{j}}$ such that $|\tilde{\theta}_{j}-\theta_{j}|\leq\epsilon$ for all $j\in[n]$ with probability $1-1/\operatorname{poly}\left(n\right)$ in time $\mathcal{O}\left(T_{U}\log{(n)}/\epsilon\right)$, where $T_{U}$ is the time to implement $U$.

Crucially, we do not need to be given access to the actual eigenvector $|v_{j}\rangle$, as this mapping can be applied to a superposition of the eigenvectors. Indeed, any quantum state $\ket*{u}$ can be decomposed in an arbitrary orthonormal basis, such as for instance in the operator eigenbasis $\ket*{v_{j}}$:

[TABLE]

where $\alpha_{j}=\bra*{u}\ket*{v_{j}}$. Hence the quantum phase estimation procedure can be applied to an arbitrary state $\ket*{u}$, which, as it is just a matter of representation, be directly applied to the operator eigenbasis, apply this mapping in the specific (eigen-)basis without knowing the actual basis:

[TABLE]

Next, we consider the process of phase kickback, where we query an oracle and encode some information pertaining to it as a relative phase.

2.6 Phase kickback

In the following, we consider a Boolean function $f:\{0,1\}^{n}\to\{0,1\}^{m}$. The function is queried via an oracle $\mathcal{O}_{f}$, and so

[TABLE]

where $|x\rangle$ is an input state and where $|q\rangle$ is an ancillary register. This operation can be implemented by a unitary circuit $U_{f}$.

If the ancillary qubit is in the state $|-\rangle=\frac{|0\rangle-|1\rangle}{\sqrt{2}}$, then by applying the oracle we obtain the state:

[TABLE]

where $\overline{\cdot}$ denotes the negated bits. By considering the cases where the output $f(x)$ is either a $1$ or a [math], it can easily be shown that this is equivalent to the state:

[TABLE]

Thus, inputs which evaluate to $1$ acquire a relative phase. This process is referred to as phase kickback. Finally, we note that we have assumed that the function can be classically computed in polynomial time, i.e. this is not an expensive step and no complexity is hidden in this call.

2.7 Amplitude amplification

Amplitude amplification [Bra+02] is an extension of Grover’s search algorithm [GR02]. Here, we shall present Grover’s search algorithm, then show how this leads to amplitude amplification.

We are given a set containing $N$ elements $\{1,\ldots,N\}$ and the goal is to find a particular element of the set which is marked. This may be modeled by the function $f:\{1,\ldots,N\}\to\{0,1\}$ such that there exists uniquely one item $a\in\{1,\ldots,N\}$ satisfying $f(a)=1$. Otherwise we have that $f(x)=0$ for all $x\in\{1,\ldots,N\}\setminus\{a\}$. Furthermore, we can evaluate $f(x)$ using an oracle that can be queried in superposition with unit computational cost, see section 2.6. Grover’s algorithm is for our purposes more appropriately seen as an algorithm for finding the (unique) root of the Boolean function $\neg f(x)$, i.e. replace our set $\{1,\ldots,N\}$ with the bitstring representations of $\{0,\ldots,N-1\}$. Then $a$ is the bitstring such that $\neg f(a)=0$, or the root of $\neg f(x)$. We shall show later how Grover’s algorithm generalises to finding all roots of a Boolean function $f(x)$.

Classically, finding $a$ requires $\mathcal{O}(N)$ oracle queries, where in the worst case we evaluate $f$ on all elements in the set. In contrast, Grover’s algorithm achieves a quadratic speed-up requiring $\mathcal{O}(\sqrt{N})$ oracle queries. Furthermore, it has been shown that there is a lower bound $\Omega(\sqrt{N})$ on the number of queries, that is, this is the optimal scaling [DH01].

Grover’s algorithm uses the phase kickback technique previously discussed in section 2.6, which associates a relative phase with the marked item. In order to do this, $n$ qubits are prepared in the $|0\rangle$ state, and by the application of the operator $H^{\otimes n}$, the uniform superposition state $|\psi_{0}\rangle=\frac{1}{\sqrt{N}}\sum_{x}|x\rangle$ is obtained. An ancillary qubit is prepared in the $|-\rangle$ state. By applying the phase kickback protocol, from Eq. (55) and (56), we obtain $\frac{1}{\sqrt{N}}\sum_{x}(-1)^{f(x)}|x\rangle|-\rangle$. We now discard the second register. The basis state corresponding to the marked item has now acquired a relative phase, and we are interested in the remaining state $|\psi\rangle=\frac{1}{\sqrt{N}}\sum_{x}(-1)^{f(x)}|x\rangle$. We denote by $Z_{f}$ the unitary operator that queries the oracle and applies the phase $(-1)^{f(x)}$ to each computational basis state $\ket*{x}$.

Ultimately, a measurement in the computational basis is to be performed. Ideally, this would yield the computational basis state corresponding to the marked item with high probability. How can this be achieved? The idea is to apply a unitary operator to the input which will dampen the coefficients associated with unmarked items, and strengthen the coefficients corresponding to the marked item. For this to be practically achievable, we need to check that this operator has an efficient implementation.

The Grover operator, to be defined shortly, does precisely this in an efficient manner. The state $|\psi\rangle$ can be expressed as a linear combination of $2^{n}$ basis terms $|x\rangle$, each with a corresponding complex coefficient $\alpha_{x}=(-1)^{f(x)}$. One of these corresponds to the marked item $a$. The idea is to strengthen the coefficient $\alpha_{a}$ whilst simultaneously weakening coefficients $\alpha_{k}$, where $1\leq k\leq N$ and $k\neq a$.

The unitary operator which achieves this is called the Grover operator and is given by $G:=DZ_{f}$, where $D=2|\psi_{0}\rangle\langle\psi_{0}|-\mathbb{I}$ is called the diffusion operator. It is not a priori obvious what this operator does or why one would wish to apply it. There are two arguments which can shed insight into this, one algebraic and one geometric, which we next introduce.

Algebraic argument. First, we consider an arbitrary state $|\psi\rangle=\sum_{x}\alpha_{x}|x\rangle$. It is easy to see that applying the diffusion operator $D$ results in the state $\sum_{x}(2\langle\alpha\rangle-\alpha_{x})|x\rangle$, where $\langle\alpha\rangle=\sum_{k}\frac{\alpha_{k}}{N}$ is defined as the mean value of the coefficients. Thus the new coefficient corresponding to $\ket*{x}$ is given by $\langle\alpha\rangle+(\langle\alpha\rangle-\alpha_{x})$. Before the first application of $D$, the mean is given by $\langle\alpha\rangle=\frac{(N-2)/N}{\sqrt{N}}$. Here, we have that the positive coefficients will be dampened and approach to zero, whereas the negative coefficients will be magnified and become positive. This is inversion about the mean. Next, the item is marked with a negative sign by $Z_{f}$ in order for the inversion about the mean to be applied in the next step. As this process repeats, the unmarked items’ coefficients will tend to zero, whereas the marked coefficient goes towards one [Wha09].

Geometric argument. Strengthening means bringing the initial state closer to the state $|a\rangle$, whilst preserving the norm. Visually, this can be seen as performing a rotation of angle $2\theta$ towards $|a\rangle$ in the plane defined by $\operatorname{span}\{\ket*{\psi_{0}},\ket*{a}\}$, where $\theta=\arcsin(\bra*{\phi_{0}}\ket*{a})$. Let $|a^{\perp}\rangle$ be the orthogonal complement of the marked item state in this subspace. Then $\theta$ is the angle between the equal superposition state $\ket*{\psi_{0}}$ and the orthogonal complement $\ket*{a^{\perp}}$.

A rotation by an angle $2\theta$ can be implemented via two successive reflections: one through $|a^{\perp}\rangle$ and one through $|\psi_{0}\rangle$. A reflection about an axis mean that any component of a particular vector orthogonal to the axis acquire a negative phase, and the component along the axis remains invariant. Let $|r\rangle$ be the axis of reflection, and let $|\psi\rangle$ be the state to reflect. After reflection, we should have $|\psi^{\prime}\rangle=\alpha_{r}|r\rangle-\sum_{i\neq r}\alpha_{i}|i\rangle$. This can be expressed as $-|\psi\rangle+2\alpha_{r}|r\rangle$, which corresponds to application of the operator $\mathbb{I}-2|r\rangle\langle r|$, up to a global phase. Thus, this gives us the reflection operation. In the case when we reflect about the axis $|\psi_{0}\rangle$ we recover the diffusion operator. In the case where we reflect about $\ket*{a^{\perp}}$ we recover $Z_{f}$. Thus, in order to implement the rotation by an angle $2\theta$, we first apply a reflection about $|a^{\perp}\rangle$, followed by a reflection about $|\psi_{0}\rangle$. We see an illustration of this in Figure 3.

We wish to have a number of rotations $k^{*}$ such that a measurement in the computational basis yields $\ket*{a}$ with high probability. This probability, after $k$ rotations, is given by $\absolutevalue{\bra*{a}G^{k}\ket*{\psi_{0}}}^{2}=\sin^{2}((2k+1)\theta)$. To have this close to one, we need $(2k^{*}+1)\theta\approx\pi/2$, which we have when $k^{*}=\frac{\pi/2}{2\theta}$. We have that $\langle\psi_{0}|a\rangle=\sin\theta$ but also $\langle\psi_{0}|a\rangle=\frac{1}{\sqrt{N}}$. Using the small angle approximation, we have that $\sin\theta\approx\theta$ and so $\theta\approx\frac{1}{\sqrt{N}}$. So, $k^{*}$ and thus Grover’s algorithm scales as $\mathcal{O}(\sqrt{N})$.

This idea can be generalized to the case where we there are now $M$ marked items out of a total of $N$ items, where $1\leq M\leq N$. First, we construct the states $|\alpha\rangle=\frac{1}{\sqrt{M}}\sum_{x_{m}}|x_{m}\rangle$ and the state over unmarked items $|\beta\rangle=\frac{1}{\sqrt{N-M}}|\sum_{x_{u}}|x_{u}\rangle$, where $\{x_{m}\}$ are the marked elements and $\{x_{u}\}$ are the unmarked elements. Once again, the goal will be to obtain a marked state with high probability when a computational basis measurement is performed. That is, the coefficients corresponding to the marked item will be boosted whereas those corresponding to unmarked items dampened. Here, we can think of this as splitting the total space $\mathcal{H}$ into two subspaces, the good subspace and the bad subspace. The good subspace corresponds to the marked basis states, i.e. those containing relevant information, which we would wish to obtain upon performing a measurement, and thus boost their corresponding amplitudes. We have that $\mathcal{H}=\mathcal{H}_{g}\oplus\mathcal{H}_{b}$, and thus we can write any arbitrary quantum state as

[TABLE]

where $|\phi_{g}\rangle$ and $|\phi_{b}\rangle$ are two orthogonal states, and where we have that $|\alpha_{g}|^{2}+|\alpha_{b}|^{2}=1$. In the context of Grover search, we have $|\phi_{g}\rangle=\frac{1}{\sqrt{M}}\sum_{x_{m}}|x_{m}\rangle$, and $|\phi_{b}\rangle=\frac{1}{\sqrt{N-M}}\sum_{x_{u}}|x_{u}\rangle$, and thus we can write the complete state as

[TABLE]

i.e. $\alpha_{g}=\sqrt{\frac{M}{N}}$ and $\alpha_{b}=\sqrt{\frac{N-M}{N}}$. The goal is to now, via the application of a unitary operator, amplify the coefficient $\alpha_{g}$ whilst weakening $\alpha_{b}$. Geometrically, this means that the quantum state will be rotated towards the good subspace.

Amplitude amplification is the process of applying Grover’s algorithm to tasks where we have an oracle for a function $f$ and need to sample $x$ from the ‘good’ subset of strings, i.e. $G=\{x\mid f(x)=1\}$, more efficiently than with a classical algorithm. By querying the oracle $\mathcal{O}\left(\sqrt{\frac{N}{M}}\right)$ times the amplitude of the ‘good’ subset of strings is then close to unity and the ‘bad’ amplitude is close to zero. Stated more generally, this gives us

Lemma 3**.**

(Amplitude Amplification)* Suppose we have an algorithm $\mathcal{A}$ that succeeds with probability $\epsilon$. Using amplitude amplification, we can take $\mathcal{O}\left(1/\sqrt{\epsilon}\right)$ repetitions of $\mathcal{A}$ to yield an algorithm $\mathcal{A}^{\prime}$ with success probability arbitrarily close to one.*

Note that $\mathcal{A}$ can be a classical or quantum algorithm and that we obtain the lemma by taking the ‘good’ subspace as successful bitstring outputs of $\mathcal{A}$ and the ‘bad’ subspace as the unsuccessful outputs. One can do this by simply appending a bit at the end of the evaluation of the result of $\mathcal{A}$ with ‘1’ for success and ‘0’ for failure.

2.8 The uncompute trick

The uncompute trick is a commonly used technique in quantum algorithms for carrying out a computation, then retrieving the initial state. For this subsection, we take much of the presentation from the discussion at [Unc]. Previously, in section 2.6, we saw that given a boolean function $f:\{0,1\}^{n}\to\{0,1\}^{m}$, there exists an oracle $\mathcal{O}_{f}$ acting as $\ket*{x}\ket*{0}\mapsto\ket*{x}\ket*{f(x)}$, which allows for the mapping $|x\rangle\to|f(x)\rangle$ to be enacted. Yet, we did not discuss how this unitary operator could be efficiently implemented. In particular, from the no-deleting theorem [KB00] there is no single-qubit unitary operator that sets an arbitrary qubit state to $\ket*{0}$. Indeed, for most algorithms we have the map $\ket*{x}\ket*{0}\ket*{0}\mapsto\ket*{x}\ket*{g(x)}\ket*{f(x)}$, where $\ket*{g(x)}$ is a garbage state in a working register we have used along the way. Here, we wish to return the garbage state to $\ket*{0}$ in the working register, so that it doesn’t disrupt future computations after we discard it. More precisely, we assume in general for any quantum algorithm that the working register is initialised in state $\ket*{0}$. If a previous computation has left the working register state as $\ket*{g(x)}$ the next algorithm to use the working register will output incorrect results. Even worse, in the case that the mapping $\ket*{x}\ket*{0}\ket*{0}\mapsto\ket*{x}\ket*{g(x)}\ket*{f(x)}$ isn’t perfect, the working register and output register are entangled. That is, we have some state $\ket*{x}\ket*{0}\ket*{0}\mapsto\sum_{y\in\{0,1\}^{m}}\alpha_{y}\ket*{x}\ket*{g(y)}\ket*{f(y)}$ and operations on the garbage register will affect the output register. Moreover, we wish to keep the input vector $\ket*{x}$. Let us see how this is achieved.

Suppose we have a unitary operation $U_{f}$ that takes an input $\ket*{x}\ket*{0}\ket*{0}$ and produces as output the state $\ket*{\phi}=\sum_{y}\alpha_{y}\ket*{x}\ket*{y}\ket*{f_{y}}$, where $\ket*{x}$ is in the input register, $\ket*{y}$ is in the working space and $\ket*{f_{y}}$ is in the output register, and $f_{y}\in\{0,1\}^{m}$, the output of the computation. Ideally, we would have that $f_{y}=f(x)$, and in this case the output is given by the state $\sum_{y}\alpha_{y}\ket*{x}\ket*{y}\ket*{f(x)}$, i.e. $\ket*{g(x)}=\sum_{y}\alpha_{y}\ket*{y}$. If an additional computational register is added, and apply the CNOT gate controlled on the third register, giving $\sum_{y}\alpha_{y}\ket*{x}\ket*{y}\ket*{f(x)}\ket*{f(x)}$. Then, if the inverse operator $U_{f}^{-1}$ is applied to the first three registers, then the state evolves to $\ket*{x}\ket*{0}\ket*{0}\ket*{f(x)}$. Finally, applying a SWAP operator111The SWAP operator, as the name suggests, swaps the state between two registers on the same number of qubits, i.e., $\operatorname{SWAP}:\ket*{x}\ket*{y}\mapsto\ket*{y}\ket*{x}$. on the last two registers produces the state $\ket*{x}\ket*{0}\ket*{f(x)}\ket*{0}$. We can now safely discard the working register as well as the ancillary register, leaving us $\ket*{x}\ket*{f(x)}$ as desired.

This action of applying $U_{f}$, appending an additional register, applying the CNOT, then applying $U_{f}^{-1}$ followed by a SWAP is known as the uncompute trick. It allows one to record the final state of a quantum computation and then reuse the working register for some other task.

By assuming $f_{y}=f(x)$ in the discussion above, we assumed that the action of $U_{f}$ was error-free, which is a highly unrealistic scenario. We now assume an error probability of $\epsilon$, that is,

[TABLE]

The output state after the uncompute operation, that is – applying the oracle, appending the ancillary register and applying the CNOT – is given by $\ket*{\phi^{\prime}}=\sum_{y}\alpha_{y}\ket*{x}\ket*{y}\ket*{f_{y}}\ket*{f_{y}}$, which is different to the ideal final state $\ket*{\phi^{\prime}}=\sum_{y}\alpha_{y}\ket*{x}\ket*{y}\ket*{f(x)}\ket*{f(x)}$. The inner product gives

[TABLE]

Finally, the inverse operator $U_{f}^{-1}$ is applied, which preserves the inner product, and thus we have that $\left\langle U_{f}^{-1}\ket*{\phi^{\prime}},U_{f}^{-1}\ket*{\phi}\right\rangle=1-\epsilon$. Thus, we have that if the unitary $U_{f}$ acts to within error $\epsilon$, the mapping $\ket*{x}\ket*{0}\mapsto\ket*{x}\ket*{f(x)}$ can be enacted to within error $\epsilon$.

2.9 Quantum RAM

In the HHL algorithm, a real-valued vector has to be manipulated by the quantum computer. That is, it must be loaded into a computational register, where the elements of the vector will be encoded in the amplitudes of a quantum state. As quantum states are normalised, these amplitudes will be the elements of the vector scaled by the norm of the vector.

More precisely, let $\mathbf{x}\in\mathbb{R}^{N}$ be the input vector and let $\ket*{x}:=\norm{\mathbf{x}}_{2}^{-1}\textstyle{\sum_{i=1}^{N}}x_{i}\ket*{i}$ be a quantum state on $\lceil\log N\rceil$ qubits. We would like to have an operation $\mathcal{R}$, that takes some finite precision representation of $\mathbf{x}$,which we call $\tilde{\mathbf{x}}$, and outputs the state $\ket*{x}$, along with $\norm{\mathbf{x}}_{2}$ i.e.

[TABLE]

where $\ket*{\widetilde{\norm{\mathbf{x}}_{2}}}$ is a computational basis state encoding the value of $\norm{\mathbf{x}}_{2}$ to some finite precision (denoted by the tilde on top).

The operation $\mathcal{R}$ is known as a quantum RAM (qRAM) and is non-trivial to implement. For instance, a linear system can be classically solved in polynomial time. Thus, a necessary condition to obtain a quantum exponential speed-up is for $\mathcal{R}$ to run in time at most polylogarithmic in $N$. This statement holds if we use a qRAM for any task that is polynomial-time computable classically.

Here, we present a memory structure used to implement it, taken from the PhD thesis of Prakash [Pra14] along with the state preparation procedure of Grover and Rudolph [GR02]. We shall say that this is a classical data structure with quantum access, in that, it stores classical information, but can be accessed in quantum superposition. This can be implemented with an ordinary quantum circuit.

First, we consider the following lemma:

Lemma 4**.**

(Controlled rotation)* Let $\theta\in\mathbb{R}$ and let $\tilde{\theta}$ be its $d$-bit finite precision representation. Then there is a unitary, $U_{\theta}$, that acts as*

[TABLE]

Proof.

Let

[TABLE]

where \sigma_{y}=\smallmatrixquantity(0&-i\\ i&0) is the Pauli $Y$ matrix. The matrix $\exp(-i\theta\sigma_{y})$ is given by \exp(-i\theta\sigma_{y})=\smallmatrixquantity(\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta). Finally, applying $U_{\theta}$ to $\ket*{\tilde{\theta}}\ket*{0}$ yields the result in Eq. (62). ∎

We should note that the unitary $U_{\theta}$ can be implemented in