Quantum Algorithm for Linear Regression
Guoming Wang∗

arXiv:1402.0660v6 [quant-ph] 31 Jul 2017

Joint Center for Quantum Information and Computer Science,
University of Maryland, College Park, MD 20742, USA
(Dated: August 1, 2017)
We present a quantum algorithm for fitting a linear regression model to a given data set using the
least squares approach. Different from previous algorithms which yield a quantum state encoding
the optimal parameters, our algorithm outputs these numbers in the classical form. So by running
it once, one completely determines the fitted model and then can use it to make predictions on new
data at little cost. Moreover, our algorithm works in the standard oracle model, and can handle
data sets with nonsparse design matrices. It runs in time poly(log(N ), d, κ, 1/ǫ), where N is the size
of the data set, d is the number of adjustable parameters, κ is the condition number of the design
matrix, and ǫ is the desired precision in the output. We also show that the polynomial dependence
on d and κ is necessary. Thus, our algorithm cannot be significantly improved. Furthermore, we
also give a quantum algorithm that estimates the quality of the least-squares fit (without computing
its parameters explicitly). This algorithm runs faster than the one for finding this fit, and can be
used to check whether the given data set qualifies for linear regression in the first place.

I.

INTRODUCTION

Curve fitting, also known as regression analysis in
statistics, is the process of constructing a mathematical function that has the best fit to a series of
data points according to some criterion. This procedure is widely used in many scientific fields, including physics, astronomy, chemistry, biology, medicine,
agriculture, geology, engineering, economics, etc. It
can help us to understand the relationship among
variables, to predict the unknown value of a variable from the known values of other variables, to
compress data, and to aid data visualization. In
practice, one often needs to fit a concise theoretical
model to a huge amount of experimental data, and
it is highly desirable to have an efficient algorithm
for this task.
Linear regression is one of the most common forms
of curve fitting. It assumes that the relationship between a dependent variable (or response) and one
or more explanatory variables (or predictors) is linear. So it fits a function which is linear in some
adjustable parameters to the given data set. These
parameters are usually determined using the (ordinary) least squares approach, which minimizes the
sum of the squared deviations of the data from the
model function. This optimization problem turns
out to be closely related to a matrix inversion problem, which is time-consuming for large data sets.
With the rise of quantum computation, one naturally asks whether quantum algorithms can perform

∗

wgmcreate@berkeley.edu

linear regression faster than their classical counterparts. Wiebe, Braun and Llyod (WBL) [1] first
studied this problem and answered it affirmatively.
Building upon the quantum algorithm for solving
linear systems of equations by Harrow, Hassidim and
Lloyd (HHL) [2], they developed a quantum algorithm for estimating the quality of the least-squares
fit for a given data set. Under the assumption that
there exist two fast procedures for specifying the
nonzero entries of the design matrix and for preparing a quantum state proportional to the response
vector, respectively (see Section II B for the definition of this matrix and vector), their algorithm has
complexity poly(log(N ), s, κ, 1/ǫ), where N is the
size of the data set, s and κ are the sparsity and
condition number of the design matrix, respectively,
and ǫ is the desired precision in the output. WBL
also gave an algorithm with similar complexity for
preparing a quantum state approximately proportional to the optimal parameters. Furthermore, they
proposed to use statistical sampling and quantum
state tomography to find a concise representation
for this state. WBL’s algorithms are mainly suited
for data sets whose design matrices are sparse and
well-conditioned.
Recently, Schuld, Sinayskiy and Petruccione
(SSP) [3] reapproached the problem of linear regression on a quantum computer from a machine learning perspective. Building upon HHL’s strategy for
matrix inversion and Lloyd, Mohseni and Rebentrost (LMR)’s density matrix exponentiation technique [4], they developed a quantum algorithm for
pattern recognition, in which one only needs to make
a prediction on a new data point based on a linear
regression model trained on a given data set and

out computing its parameters explicitly). This algorithm runs faster than the one for finding this fit,
and can be used to check whether the given data set
qualifies for linear regression in the first place.
We make use of two recent results in designing
our algorithms. The first one is Low and Chuang’s
method for Hamiltonian simulation based on qubitization [5] and quantum signal processing [6]. This
method allows us to simulate a nonsparse Hamiltonian, provided that this Hamiltonian can be embedded into a larger unitary operator in certain
way. The second one is Childs, Kothari and Somma
(CKS)’s approach to matrix inversion [7]. This approach differs from HHL’s in that it does not use
phase estimation, but relies on a techique for implementing a linear combination of unitaries (LCU)
and a suitable Fourier or Chebyshev series representation of the matrix inverse function. Consequently,
it has exponentially better dependence on the precision than HHL’s approach. We combine these results with traditional techniques (such as amplitude
estimation [8]) to find the optimal parameters and
to estimate the quality of the least-squares fit for a
given data set.
As mentioned before, WBL have suggested a
sampling-based algorithm for learning the optimal
parameters in Ref. [1]. Our algorithm for computing
the optimal parameters differs from their algorithm
in several ways. First, as mentioned above, our algorithm uses the approach of Ref. [7] for matrix inversion, which has better dependence on the desired
precision in the output than HHL’s approach (which
was used by Ref. [1]). Second, we compute the pseudoinverse of the design matrix by considering its
singular value decomposition (SVD), while Ref. [1]
achieved this by following a step-by-step approach
(see the end of Section IV for more discussion on
this). Third, as mentioned above, our algorithm is
based on the method of Ref. [5] for simulating a large
class of Hamiltonians, while Ref. [1] was based on an
old method for simulating sparse Hamiltonians. As
a consequence, our algorithm can handle data sets
with nonsparse design matrices. Fourth, we assume
that the data set is given via standard oracles (see
Section II C for more details), and explicitly address
the issue of preparing a quantum state proportional
to the response vector. By contrast, Ref. [1] ignored
the cost of this step. Finally, our algorithm uses amplitude estimation to estimate the optimal parameters, which has quadratically better dependence on
the desired accuracy in the output than statistical
sampling (which was used by Ref. [1]).
The remainder of this paper is organized as follows. In Section II, we provide some requisite back-

does not need to find this model explicitly. Their
algorithm takes as input multiple copies of three
quantum states encoding the design matrix of the
training set, the response vector of the training set,
and the new data point, respectively, and outputs
a scalar value which is the predicted response for
the new data point. Excluding the costs of preparing these states and assuming the design matrix is
close to a low-rank matrix, this algorithm has complexity poly(log(d), κ, 1/ǫ), where d is the number
of adjustable parameters, κ is the condition number
of the design matrix, and ǫ is the desired precision
in the output. SSP’s algorithm is mainly suited for
data sets whose design matrices are well-conditioned
and have low-rank approximations.
Both WBL and SSP have focused on the scenario
where both the size N of the data set and the number
d of adjustable parameters are exponentially large.
Thus, they do not attempt to find the optimal parameters explicitly (which is time-consuming), but
only encode these parameters in a quantum state
(which can be used to make predictions on new data
via swap test). While this scenario is useful in some
applications (e.g. estimation of the output state of a
quantum device), we believe that it is equally important to consider the scenario where d is much smaller
than N . Namely, N is exponentially large, but d is
only polynomially large. One often encounters this
situation when dealing with a classical data set and
wanting to compress a large amount of data into a
concise model (with few parameters). Once such a
model is found explicitly, one can use it to make
predictions on new data at little cost. Furthermore,
saving the optimal parameters is much easier than
storing the quantum state encoding these numbers,
as quantum resources are fragile.
For the above reasons, in this paper, we present a
new quantum algorithm for fitting a linear regression
model to a given data set using the least squares approach. Our algorithm works in the standard oracle
model, and outputs the optimal parameters in the
classical form. It runs in time poly(log(N ), d, κ, 1/ǫ),
where N is the size of the data set, d is the number
of adjustable parameters, κ is the condition number
of the design matrix, and ǫ is the desired precision
in the output. Note that the polynomial dependence
on d is inevitable, because simply writing down all
the optimal parameters takes Ω(d) time. We show
that the polynomial dependence on κ is also necessary, by proving a lower bound on the quantum
query complexity of this problem. These facts imply
that our algorithm cannot be significantly improved.
Furthermore, we also give a quantum algorithm that
estimates the quality of the least-squares fit (with2

norm than the quadratic mean norm of x’s entries.
In particular, we say that x is balanced if ρ(x) =
O(1) (e.g. at most 100).
Given a matrix A = (ai,j ) ∈ CN ×M , we define ai := (ai,1 , ai,2 , . . . , ai,M )T , for each i ∈
{1, 2, . . . , N }. We also use kAk and kAkF to denote
the spectral and Frobenius norms of A, respectively,
i.e.
kAxk
kAk :=
max
,
(4)
M
x∈C , x6=0 kxk

ground information, and formally state the problems
studied in this work. In Section III, we describe
an efficient procedure for simulating a nonsparse
Hamiltonian related to the design matrix, which is
a key component of our algorithms. In Section IV,
we present a quantum algorithm for fitting a linear regression model to a given data set using the
least squares approach. In Section V, we propose a
quantum algorithm that estimates the quality of the
least-squares fit (without computing its parameters
explicitly). In Section VI, we prove a lower bound
on the quantum query complexity of linear regression. Finally, we conclude in Section VII with some
comments and future research directions.

and

v
uN M
uX X
2
|ai,j | .
kAkF := t

(5)

i=1 j=1

II.

In addition, we define

PRELIMINARIES

kAk2,∞ :=

In this section, we provide the necessary background information to understand this paper. In
Section II A, we introduce the notation used in this
paper. In Section II B, we review some basic facts
about linear regression. In Section II C, we formally
state the problems studied in this work.

1≤i≤N

Notation

σ(A) :=

Given a real number x, we define its sign as
sgn(x) = 1 if x ≥ 0, and sgn(x) = −1 otherwise.
Given two real numbers a, b and a real number δ > 0,
we say that a is a δ-additive approximation of b if
|a − b| ≤ δ. Moreover, we say that an algorithm
estimates a quantity x up to additive error δ if it
outputs a δ-additive approximation of x.
Given a vector x = (x1 , x2 , . . . , xN )T ∈ CN , we
use kxk∞ and kxk to denote the l∞ and l2 norms of
x, respectively, i.e.
kxk∞ := max |xi |,

(1)

v
uN
uX 2
kxk := t |xi | .

(2)

1≤i≤N

√
N kAk2,∞
kAkF

(6)
(7)
(8)

j=1

max kai k
1≤i≤N
= s
.
N
P
2
1
kai k
N

(9)

i=1

The smaller σ(A) is, the more balanced A is, in the
sense that no row of A has significantly larger norm
than the quadratic mean norm of A’s rows. In particular, we say that A is balanced if σ(A) = O(1)
(e.g. at most 100).
For the above A, we also use Range(A) to denote the range (i.e. column space) of A, and use
Π(A) to denote the projection onto Range(A). We
also use sj (A) to denote j-th smallest singular value
of A (counted with multiplicity), and use λj (A) to
denote the j-th smallest eigenvalue of A (counted
with multiplicity), starting with j = 1. The condition number of A, denoted by κ(A), is defined
as the ratio of largest to smallest singular value of
A. Futhermore, we use A+ to denote the MoorePenrose pseudoinverse of A. That is, if A has
P
the singular value decomposition A = k sk uk vk† ,
where sk > 0, uk ∈ CN and vk ∈ CM are unit
P
†
vectors, then A+ := k s−1
k vk uk .
N ×M
Given a matrix A ∈ C
and a vector x ∈ CN ,
we define

and

i=1

Moreover, we define
√
max |xi |
N kxk∞
1≤i≤N
.
ρ(x) :=
= s
kxk
N
P
2
1
|xi |
N

kAxk∞
kxk

= max kai k
1≤i≤N
v
uM
uX
2
= max t |ai,j |

and
A.

max

x∈CM , x6=0

(3)

i=1

The smaller ρ(x) is, the more balanced x is, in the
sense that the no entry of x has significantly larger

τ (A, x) :=

3

kΠ(A)xk
kxk2

2

.

(10)

Vm−1 → Vm . Suppose Pi is a procedure that implements Vi with precision ǫi and failure probability
O(ǫi ), for each i ∈ {1, 2, . . . , m}. Let P be the concatenation of these procedures (i.e. P1 → P2 →
· · · → Pm−1 → Pm ). Then by a standard hybrid
argument, one can show that P implements the unitary P
operation V := Vm Vm−1 . . . V2 V1 with precision
ǫ := m
i=1 ǫi and failure probability O(ǫ). Thus, on
any input state |ψi, the procedure P succeeds with
probability 1 − O(ǫ) (with a flag indicating success),
and when it succeeds, it outputs a state O(ǫ)-close
to V |ψi in l2 norm. This fact will be useful in the
design of our algorithms.

In words, τ (A, x) measures how much “fraction” of
x lies in the range of A. In particular, we say that
(A, x) is well-behaved if τ (A, x) = Ω(1) (e.g. at least
2/3).
Given a vector x ∈ CN , we say that x is d-sparse
if it contains at most d nonzero entries. Given a
matrix A ∈ CN ×M , we say that A is d-sparse if it
contains at most d nonzero entries in each row and
column. In particular, if d = poly(log(L)) where
L = max{N, M }, then we simply say that A is
sparse.
Given a state |ϕi ∈ Cd and a real number ǫ > 0,
we say that a procedure prepares |ϕi with precision ǫ
if this procedure prepares a state |ψi ∈ Cd satisfying
k|ϕi − |ψik ≤ ǫ.
Given a unitary operation V ∈ U(d) and a real
number ǫ > 0, we say that a procedure implements
V with precision ǫ and failure probability O(ǫ) if
there exists an integer l ≥ 0 such that, on any input
state |ψi ∈ C(d), this procedure first appends an lqubit ancilla system in state |0l i, then performs a
unitary operation U ∈ U(2l × d) on the joint system
such that
X
|jiBj |ψi,
(11)
U 0l |ψi = 0l A|ψi +

B.

Linear Regression
N

Given a data set {yi , xi,1 , xi,2 , . . . , xi,d }i=1 of N
statistical units (where N ≥ d), a linear regression
model assumes that the relationship between the response (or regressand, dependent variable) yi and
the predictors (or regressors, explanatory variables)
xi,1 , xi,2 , . . . , xi,d is linear. That is, there exist some
unknown parameters β1 , β2 , . . . , βd and residual
terms ǫi such that

j6=0l

yi = β1 xi,1 + β2 xi,2 + · · · + βd xi,d + ǫi , 1 ≤ i ≤ N.
(18)
In the matrix form, it can be written as

where A and the Bj ’s are linear operators satisfying
P
kA − V k ≤ ǫ and A† A+ j6=0l Bj† Bj = I, and finally
measures the ancilla system and postselects on the
outcome being 0l . Note that since V is unitary and
kV − Ak ≤ ǫ, we get kV |ψi − A|ψik ≤ ǫ, and
|kA|ψik − 1| = |kA|ψik − kV |ψik|
≤ k(A − V )|ψik
≤ ǫ,

y = Xβ + ǫ,
where

 
xT1
x1,1 x1,2 . . . x1,d
 xT2   x2,1 x2,2 . . . x2,d 

  
X :=  .  =  .
. ,
.. . .
 ..   ..
. .. 
.
xN,1 xN,2 . . . xN,d
xTN
 
 
 
ǫ1
β1
y1
 ǫ2 
 β2 
 y2 
 
 
 
y :=  . , β :=  . , ǫ :=  . .
 .. 
 .. 
 .. 
ǫN
βd
yN


(12)
(13)
(14)

and
V |ψi −

A|ψi
≤ kV |ψi − A|ψik
kA|ψik
A|ψi
+ A|ψi −
kA|ψik
≤ ǫ + |kA|ψik − 1|
≤ 2ǫ.

(19)

(15)

(20)

(21)

We usually call X the design matrix, y the response
vector, β the parameter vector, and ǫ the residual
vector. Here we assume that the xi,j ’s and yi ’s are
real numbers. This is actually without loss of generality, because any linear regression model with complex variables can be reduced to a (slightly larger)
linear regression model with real variables. Moreover, we assume that the design matrix X has full
rank d. In other words, the d columns of X are linearly independent. This is a necessary condition for
linear regression to have a unique solution.

(16)
(17)

Thus, on any input state |ψi, this procedure suc2
ceeds with probability kA|ψik = 1 − O(ǫ) (with a
flag indicating success), and when it succeeds, it outA|ψi
which is O(ǫ)-close to V |ψi in
puts the state kA|ψik
2
l norm.
Now consider a quantum circuit consisting of a
sequence of unitary operations V1 → V2 → · · · →
4

Then we have

We emphasize that the predictors can be nonlinear
functions of some “baseline” variables. This allows
linear regression to fit a nonlinear relationship between the response and the baseline variables. For
example, suppose we are interested in learning how
the yield yi of a chemical synthesis is related to the
temperature ti at which the synthesis takes place.
We propose a quadratic model of the form:
yi = a0 + a1 ti + a2 t2i + ǫi , 1 ≤ i ≤ N.

+  
0 X
y
0
XT 0
 

T +
y
0 (X )
=
0
X+
0


0
=
X+ y
 
0
=
β̂
= z.

A+ b =

(22)

This model is linear in the parameters a0 , a1 and
a2 , but nonlinear in the baseline variable ti . In the
matrix form, it can be written as
 

  
ǫ1
1 t1 t21  
y1
 ǫ2 
 y2  1 t2 t22  a0
 

  
 ..  =  .. .. .. a1  +  ..  (23)
 . 
 .  . . .  a
2
ǫN
yN
1 tN t2N

(29)
(30)
(31)
(32)
(33)

The fact that A is a real symmetric matrix facilitates
the implementation of A+ . Once we have a procedure for preparing a quantum state proportional to
A+ b = z, we can utilize this procedure to get useful
information about β̂.
A statistical model fits a data set well only if the
discrepancy between the observed response and the
response predicted by this model is small. Here we
measure the quality of the least-squares fit y ≈ Xβ̂
using the quantity

Here the design matrix is a Vandermonde matrix,
and it has full rank as long as there are at least
three distinct ti ’s. Furthermore, this design matrix
is not sparse. This is a generic phenomenon in linear
regression, because we often include the constant 1
as one of the predictors, and consequently the design
matrix often contains a dense column of all 1’s.
Linear regression models are usually fitted using
the least squares approach, which minimizes the sum
of the squared residuals. Namely, it finds
2
β̂ := argmin kXβ − yk



τ :=

kŷk2

kyk

2

=1−

kǫ̂k2

kyk

2,

(34)

where

(24)

β∈Rd

ŷ := Xβ̂ = Π(X)y

(35)

ǫ̂ := y − ŷ = (I − Π(X))y.

(36)

and

This optimization problem has the following closedform solution [9]
β̂ = X+ y = (XT X)−1 XT y.

(25)

Π(X) = X(XT X)−1 XT ,

(26)

Xβ̂ = Π(X)y.

(27)

Namely, 1 − τ is the ratio of the squared norm of
the residual vector ǫ̂ to that of the response vector
y. It turns out that τ = τ (X, y). So a data set
(X, y) can be explained well by a linear regression
model only if it is well-behaved, i.e. τ (X, y) = Ω(1)
(i.e. at least 2/3). This kind of data sets will be the
main focus of our study. We will also give an efficient
quantum algorithm for testing whether a given data
set is well-behaved or not.
It is worth noting that Wiebe, Braun and
Lloyd (WBL) [1] have used the quantity E =
2
k|yi − F|λik (according to their notation) to measure the error of the least-squares fit. Their |yi, F
and |λi correspond to our y/kyk, X and β̂/kyk, respectively. Then one can see that their 1 − E is
equivalent to our τ . So WBL essentially measured
the quality of the least-squares fit in the same way
as we do.
In practice, after one collects the raw data from
the experiements, one does not immediately fit a

Noting that
we obtain

Namely, Xβ̂ is exactly the projection of y onto the
range of X. This is the geometric interpretation of
least-squares linear regression.
Although Eq. (25) gives the solution of linear regression, it is not computationally convenient, because X is a rectangular matrix and X+ is not easy
to implement physically. To overcome this issue, we
adopt the strategy of Ref. [2] (which was also used
in Ref. [1]) and embed β̂ into the solution of a larger
linear system. Specifically, let
 


 
0
0 X
y
. (28)
, z :=
A :=
, b :=
0
XT 0
β̂
5

there exists a procedure Px that allows us to perform the map

mathematical function to these data. Instead, one
needs to preprocess the raw data to make them wellsuited for data fitting. This preprocessing usually
consists of imputation of missing data, data normalization or standadization, and elimination of influential outliers which have detrimental effect on the
estimated regression function (an outlier is a data
point whose response y does not follow the general
trend of the rest of the data). The last step is important, because we want the fitted model to capture the typical relationship among the response and
predictors so that it can be generalized to new data.
This requires that the loss function
2

kXβ − yk =
=

N
X
i=1

N
X
i=1

xTi β − yi



d
X
j=1

2

|ii|ji|zi 7→ |ii|ji|z ⊕ xi,j i

for any i ∈ {1, 2, . . . , N } and j ∈ {1, 2, . . . , d}, where
the third register holds a bit string representing an
entry of X. For y = (y1 , y2 , . . . , yN )T ∈ RN , we
assume there exists a procedure Py that allows us to
perform the map
|ii|zi 7→ |ii|z ⊕ yi i

2

βj xi,j  − yi ,

(38)

should not be dominated by only a few data points.
An outlier has the potential to do so, especially if
it has high leverage (i.e. it has “extreme” predictor
x values). However, we emphasize that not all outliers are influential and should be eliminated. The
identification of influential outliers is an important
and complicated topic, and many techniques have
been developed for this task, such as difference in
fits (DFFITS) and Cook’s distance. In this paper,
we assume that the given data set has already been
preprocessed, and the harmful outliers have been removed. Since there is no general characterization
of such data points, we assume for simplicity that
no xi or yi has extremely large norm (compared to
the average norm of the xi ’s or yi ’s, respectively).
In other words, both X and y are balanced, i.e.
σ(X) = O(1) (e.g. at most 100) and ρ(y) = O(1)
(e.g. at most 100). These assumptions ensure that
no data point has significantly larger contribution
to the loss function than the others, and are useful in practice. However, we acknowledge that these
assumptions might be too stringent for some applications, and standard preprocessing techniques do
not always guarantee them, and it remains future
work to extend our results to the most general case
of linear regression.
C.

(40)

for any i ∈ {1, 2, . . . , N }, where the second register
holds a bit string representing an entry of y. We assume that both Px and Py are efficient, in the sense
that they run in time poly(log(N )). This requires
that either each entry of X and y can be quickly
computed by an algorithm (given its position), or
these entries are stored in a quantum random access
memory (QRAM) beforehand. Our algorithms work
well in both cases.
Given access to Px and Py , our primary goal
is to fit a linear regression model to the data set
{yi , xi,1 , xi,2 , . . . , xi,d }N
i=1 using the least squares approach. Our secondary goal is to estimate the quality of the fitted model (without computing its parameters explicitly).
Formally, we define our linear regression (LR)
problems as follows:

(37)


(39)

Problem 1 (LR-P). Let X = (xi,j ) ∈ RN ×d be a
balanced matrix such that its singular values are in
the range [1/κ, 1]. Let y = (y1 , y2 , . . . , yN )T ∈ RN
be a balanced unit vector. Suppose (X, y) is wellbehaved. Given ǫ > 0 and access to the procedures Px
and Py described above, the goal is to output a vector
β := (β1 , β2 , . . . , βd )T ∈ Rd satisfying kβ − β̂k∞ ≤ ǫ,
where β̂ := X+ y, succeeding with high probability
(e.g. at least 2/3).
Problem 2 (LR-Q). Let X = (xi,j ) ∈ RN ×d be a
balanced matrix such that its singular values are in
the range [1/κ, 1]. Let y = (y1 , y2 , . . . , yN )T ∈ RN
be a balanced unit vector. Given ǫ > 0 and access to the procedures Px and Py described above,
the goal is output an ǫ-additive approximation of
τ := kΠ(X)yk2 /kyk2 , succeeding with high probability (e.g. at least 2/3).

Problem Statement

Although in the above problems we assume that
the singular values of X lie in the range [1/κ, 1] and
kyk = 1, this is without loss of generality. Suppose
instead that the singular values of X lie in the range
[a/κ, a] and kyk = b, for some constants a, b > 0.

In this paper, we assume that the data set
N
{yi , xi,1 , xi,2 , . . . , xi,d }i=1 is given via two black-box
subroutines. For X = (xi,j ) ∈ RN ×d , we assume
6

Namely, X and y are rescaled by a factor of a and
b, respectively. Then β̂ = X+ y is rescaled by a
factor of b/a. So we only need to multiply the result
of LR-P by this factor. On the other hand, τ =
kΠ(X)yk2 /kyk2 is immune to this rescaling. So we
do not need to make any change to the result of
LR-Q.
We will develop quantum algorithms for solving
the above problems. We quantify the resource requirements of these algorithms using two measures.
The query complexity is the number of uses of the
procedures Px and Py in the algorithm. The gate
complexity is the number of 2-qubit gates used in
the algorithm. An algorithm is gate-efficient is if
it is gate complexity is larger than its query complexity only by a logarithmic factor. Formally, an
algorithm with query complexity Q is gate-efficient
if its gate complexity is O(Q · poly(log(QN ))). All
the algorithms presented in this paper will be gateefficient.
III.

Lemma 1. Let X be defined as in LR-P or LR-Q.
Let A := |1ih0| ⊗ XT + |0ih1| ⊗ X. Then there exists
a gate-efficient procedure that simulates e−iAt with
precision ǫ and failure probability O(ǫ) by making
 
 
√
1
O d
dt + log
ǫ
uses of Px .

√
√
Proof. Let X̃ = (x̃i,j ) := X · N /(σ d), where σ :=
σ(X) = O(1). Then we claim
X̃

2,∞

= max kx̃i k
1≤i≤N
v
u d
uX
2
= max t |x̃i,j |
1≤i≤N

(41)
(42)

j=1

≤ 1.

(43)

To see this, recall that the singular values of X are
in the range [1/κ, 1]. So

HAMILTONIAN SIMULATION
2
kXkF

Hamiltonian simulation is an important topic that
has received a lot of attention in the past years
[5, 6, 10–21]. Recently, Low and Chuang [6] proposed a techinque named quantum signal processing,
and showed how to use this technique and Childs’
quantum walk [14] to simulate sparse Hamiltonians
nearly optimally. Later, in Ref. [5], they proposed
another technique called qubitization, and demonstrated how to use this technique and quantum signal processing to simulate a larger class of Hamiltonians efficiently. In this paper, we will use the
following variant of their result:

d
 X
(sj (X))2 ≤ d.
= tr X X =
T

(44)

j=1

This implies that
kXk2,∞ = max kxi k
1≤i≤N

(45)

σkXkF
√
(46)
N
√
σ d
(47)
≤ √ .
N
√
√
Using this fact and X̃ = X · N /(σ d), we obtain
Eq. (43), as desired.
Now let V̂ be a unitary operator such that
=

Theorem 1 (Adapted from Theorem 1 of Ref. [5]).
Let Û and Ĝ be unitary operators on n and k (< n)
qubits, respectively, such that hG|Û |Gi = Ĥ is a Hermitian operator on n−k qubits, where |Gi := Ĝ 0k .
Then there exists a gate-efficient algorithm that simulates e−iĤt with precision ǫ and failure probability
O(ǫ) by making O(t + log(1/ǫ)) uses of controlled-Ĝ
and controlled-Û .

V̂ |0, ii1 |0, 0m i2 |0i3 = |0, ii1 |ϕi i2,3 , 1 ≤ i ≤ N,
V̂ |1, ji1 |0, 0m i2 |0i3 = |1, ji1 |ψi2,3 , 1 ≤ j ≤ d,

(48)

(49)

where m = Θ(log(N )),

Theorem 1 provides a way to simulate a nonsparse
Hamiltonian, provided that this Hamiltonian can be
embedded into a larger unitary operator in the way
described above. Using this fact, we develop an efficient procedure for simulating e−iAt , which will be a
crucial component of our algorithms for solving the
LR-P and LR-Q problems. Recall that A is defined
by Eq. (28) and is not sparse in general.

|ϕi i2,3

:=

d
X
j=1

q
2
x̃i,j |1, ji2 |0i3 + 1 − kx̃i k |1, 1i2 |1i3

(50)

and
N
1 X
|ψi2,3 := √
|0, ii2 |0i3 .
N i=1

7

(51)

O(N )-dimensional Hilbert space [22]. Finally, we
uncompute the xi,j ’s in the last register by making O(d) uses of Px , and obtain the desired state
|0, ii1 |ϕi i2,3 . This process requires O(d) uses of Px
and is gate-efficient. Combining the above facts,
we know that V̂ can be implemented by a gateefficient procedure that makes O(d) uses of Px , as
claimed.

Let SWAP1,2 be the swap operator on the first two
registers, i.e. SWAP1,2 |ϕi1 |ψi2 = |ψi1 |ϕi2 for all
states |ϕi and |ψi. Then we define
Ŵ := (SWAP1,2 ⊗I3 ) · V̂

(52)

Û := Ŵ † V̂ .

(53)

and
In addition, let |Gi := |0, 0m i2 |0i3 . Then by a direct
calculation, one can verify that
A
Ĥ := hG|Û |Gi = √ .
σ d

We remark that our embedding construction in
the proof of Lemma 1 looks similar to the construction in Refs. [14, 18]. However, we emphasize that
our high-level strategy for simulating the Hamiltonian A is very different from that of Refs. [14, 18].
Specifically, Refs. [14, 18] simulate a Hamiltonian
by embedding it into a quantum walk operator and
then performing phase estimation on this operator.
By contrast, we simulate e−iAt by embedding A
into a unitary operator Û (which is not a quantum
walk) in certain way and then invoking the method
of Ref. [5] for Hamiltonian simulation (which is arguably more advanced than phase-estimation-based
methods). So the similarity between the proof of
Lemma 1 and the construction in Refs. [14, 18] is
superficial rather than essential.

(54)

We will show below that Û can be implemented by
a gate-efficient procedure that makes use O(d) uses
of Px . Then by Theorem 1,
e−iAt = e−iĤσ

√

dt

(55)

can be implemented with precision ǫ and failure
probability O(ǫ) by a gate-efficient procedure that
makes
 
 
√
1
O d σ dt + log
(56)
ǫ
 
 
√
1
=O d
dt + log
(57)
ǫ

IV.

uses of Px (recall that σ = O(1)), as claimed.
Clearly, SWAP1,2 can be implemented in time
poly(log(N )). So it remains to show that V̂ can
be implemented by a gate-efficient procedure that
makes use O(d) uses of Px . To prove this, first note
that the mapping
|1, ji1 |0, 0m i2 |0i3 → |1, ji1 |ψi2,3

In this section, we present a quantum algorithm
for solving the LR-P problem, i.e. finding the parameters β̂ = X+ y of the least-squares fit y ≈ Xβ̂
for a given data set (X, y). Roughly speaking, this
algorithm computes β̂ = (β̂1 , β̂2 , . . . , β̂d )T in three
stages. The first stage estimates the absolute values
of the β̂i ’s. The second stage determines the signs
of these parameters, up to a global sign ±1. That
is, up to this stage, we obtain a vector β ∈ Rd which
is close to either β̂ or −β̂. The final stage decides
which of the two cases holds. This algorithm relies
on several subroutines (besides the one for Hamiltonian simulation in Lemma 1). One is the following
PN
procedure for preparing the state |yi = i=1 yi |ii
(recall that kyk = 1).

(58)

can be implemented in time O(log(N )), since
PN
|ψi2,3 = √1N i=1 |0, ii2 |0i3 is easy to prepare.
Meanwhile, we can accomplish the transformation
|0, ii1 |0, 0m i2 |0i3 → |0, ii1 |ϕi i2,3

(59)

as follows. First, we learn xi,1 , xi,2 , . . . , xi,d by
making O(d) uses of Px , and obtain the state


d
O
|0, ii |0, 0m i |0i 
|xi,j i .
(60)
1

2

3

j=1

FINDING THE LEAST-SQUARES FIT

Lemma 2. Let y be defined as in LR-P or LR-Q.
P
Then the state |yi = N
i=1 yi |ii can be prepared with
precision δ by a gate-efficient procedure that makes
O(log(1/δ)) uses of Py .

4

Then, we perform a unitary operation on the second and third registers depending on the content of the last register, and convert |0, 0m i2 |0i3
into |ϕi i2,3 . This step can be achieved in time
Pd
O(d · log(N )), since |ϕi i2,3 = j=1 x̃i,j |1, ji2 |0i3 +
q
1 − kx̃i k2 |1, 1i2 |1i3 is a O(d)-sparse vector in an

Proof. Consider the following procedure which
transforms |0n i into |yi probabilistically, where n =
PN
Θ(log(N )). First, we map |0n i to √1N i=1 |ii in
time O(log(N )). Then, we convert this state into
8

√1
N

PN

Proof. Let A := |1ih0| ⊗ XT + |0ih1| ⊗ X and |bi :=
|0i|yi. Suppose X has the singular value decomposition

i=1 |ii|yi i

by making O(1) uses of Py . Next,
we append an ancilla qubit in state |0i, and perform
the controlled-rotation


s
2
y
|y
|
i
i
|yi i|0i → |yi i
|0i + 1 −
|1i (61)
kyk∞
kyk2∞

X=

j=1

A=

l

n

l

A 0 |0 i − 0 |yi

(65)

sj |+j ih+j | −

d
X
j=1

sj |−j ih−j |,

(66)

where
1
|±j i := √ (|0i|uj i ± |1i|vj i).
2

(67)

So A is a Hermitian matrix whose nonzero eigenvalues are in the range Dκ := [−1, −1/κ] ∪ [1/κ, 1].
Moreover, |bi = |0i|yi is a unit vector, and A+ |bi =
|1i|β̂i by Eq. (33).
We will use a recent technique proposed by Childs,
Kothari and Somma [7] to approximately invert the
matrix A. Let the function h(x) be defined as

√
√
1 − δ ′ 0l |yi + δ ′ Φ⊥ , (62)

2

d
X
j=1

h(x) :=


where l is a positive integer, δ ′ = O δ 2 , and Φ⊥ is
a normalized state satisfying ( 0l ih0l ⊗ I) Φ⊥ = 0.
This implies that
′

sj |uj ihvj |,

where sj ∈ [1/κ, 1], |uj i ∈ RN and |vj i ∈ Rd are
unit vectors, for all j ∈ {1, 2, . . . , d}. Then A has
the spectral decomposition

onthe last two registers, where kyk∞ = maxi |yi | =
Θ √1N (since y is a balanced unit vector). After
that, we measure the ancilla qubit, and with proba2
bility Nkyk
kyk2∞ = Ω(1), the outcome is 0 and we obtain
PN
the state i=1 yi |ii|yi i. Finally, we uncompute yi in
the second register by making O(1) uses of Py , and
PN
obtain the desired state |yi = i=1 yi |ii.
The above procedure, denoted by A, makes O(1)
uses of Py , is gate-efficient, and has success probability Ω(1). We use Grover’s π/3-amplitude amplification (i.e. the generalization of fixed-point quantum search)
 [23] to raise the success probability to
1 − O δ 2 . This boosted procedure, denoted by A′ ,
requires O(log(1/δ)) repetitions of A, and satisfies
A′ 0l |0n i =

d
X

J−1
X

K
X

α(j, k)e−ixη(j,k) ,

(68)

j=0 k=−K

where
i
2 2
α(j, k) := √ kδy δz2 e−k δz /2 ,
2π
η(j, k) := jkδy δz ,

√
= (1 − 1 − δ ′ )2 + δ ′ (63)

(64)
= O δ2 .

(69)
(70)

for some J
= Θ((κ/ǫ)
K
=
 p· log(κ/ǫ)),

Θ(κ · log(κ/ǫ)), δy = Θ ǫ/ log(κ/ǫ) and δz =


p
Θ 1/(κ log(κ/ǫ)) . Then h(x) is ǫ-close to 1/x

Furthermore, A′ makes O(log(1/δ)) uses of Py , and
is gate-efficient. So A′ satisfies all the desired properties. This concludes the proof.

on the domain Dκ [7], i.e.

Our algorithm for solving the LR-P problem also
requires the following procedures for computing |β̂i |
and |β̂i − β̂j |.

h(x) − x−1 ≤ ǫ,

∀x ∈ Dκ .

(71)

Then since A is a Hermitian matrix whose nonzero
eigenvalues are in the range Dκ , we have

Lemma 3. Let X, y and β̂ be defined as in LR-P.
Then there exists a gate-efficient quantum algorithm
that makes
 1.5 3
  κ 
d κ
·
poly
log
O
ǫ2
ǫδ

h(A) − A+ ≤ ǫ.

(72)

This implies that
h(A)|bi − A+ |bi ≤ ǫ,

(73)

as |bi is a unit vector. Moreover, Ref. [7] shows that

uses of Px and Py , and outputs an ǫ-additive approximation of |β̂i |, for any given i ∈ {1, 2, . . . , d},
succeeding with probability at least 1 − δ.

α :=

J−1
X

K
X

j=0 k=−K

9

 p

|α(j, k)| = Θ κ log(κ/ǫ) ,

(74)

+ Ξ⊥

and
|η(j, k)| ≤ JKδy δz = Θ(κ · log(κ/ǫ)),

(75)

≤ k|z′ i − |zik
= O(ǫ).

p′ :=

(77)
(78)

(79)

where m = O(log(JK)) = O(log(κ/ǫ)), and let U
be defined as
J−1
X

K
X

j=0 k=−K

|j, kihj, k| ⊗ sgn(k)e−iAη(j,k) . (80)

(81)

A direct calculation shows that
1
W |0m i|bi = |0m ih(A)|bi + Φ⊥
(82)
α


kh(A)|bik
h(A)|bi
=
|0m i
α
kh(A)|bik
+ Φ⊥ ,

(83)

where Φ⊥ is an unnormalized state satisfying
(|0m ih0m | ⊗ I) Φ⊥ = 0. Next, let R be a unitary
operator such that

repetitions of R, W = V † U V and the procedure
for preparing |bi = |0i|yi. This means that we
need to implement U with precision O(ǫ′′ ) and failure probability O(ǫ′′ ). We also need to prepare |yi
with precision O(ǫ′′ ). By Lemma 1, Eqs. (75) and
(80), and Lemma 8 of Ref. [7], U can be implemented with precision O(ǫ′′ ) and failure probability O(ǫ′′ ) by a gate-efficient
procedure that makes

O d1.5 κ · poly(log(κ/ǫ)) uses of Px . Meanwhile, by
Lemma 2, |yi can be prepared with precision O(ǫ′′ )
by a gate-efficient procedure that makes O(log(κ/ǫ))
uses of Py . Furthermore, V can be implemented in
time O(κ · poly(log(κ/ǫ))) [7], and clearly R can be

R|0i|1, ii = |1i|1, ii,
(84)
R|0i|1, i′ i = |0i|1, i′ i, 1 ≤ i′ ≤ N, i′ 6= i, (85)
R|0i|0, ji = |0i|0, ji, 1 ≤ j ≤ N.
(86)
Then by Eqs. (83), (84), (85) and (86), we obtain
h1, i|z′ i m
|0 i1 |1i2 |1, ii3
α
X h1, i′ |z′ i
+
|0m i1 |0i2 |1, i′ i3
α
′

RW |0m i1 |0i2 |bi3 =

i 6=i

+

X h0, j|z′ i
j

α

(88)

= O(ǫ).
(92)
√ ′
Namely, α p̂ is an O(ǫ)-additive approximation of
|h1, i|z′ i|, as desired.
The above basic algorithm has success probability
at least 3/4. To boost the success probability to at
least 1 − δ, we repeat this algorithm O(log(1/δ))
times, and take the median of the estimates from
these runs. A standard Chernoff’s bound ensures
that the failure probability is at most δ.
Let us analyze the complexity of this algorithm.
Since we want to estimate p′ up to additive error ǫ′′ ,
amplitude estimation requires
 
 2
 2

1
α
κ log(κ/ǫ)
O ′′ = O 2 = O
(93)
ǫ
ǫ
ǫ2

Then we define
W := V † U V.

|h1, i|z′ i|
.
α2

We use amplitude estimation [8] to obtain an ǫ′′ additive approximation p̂′ of p′ , where
 2


ǫ
ǫ2
′′
ǫ := Θ 2 = Θ 2
,
(89)
α
κ log(κ/ǫ)
√
succeeding
Then p̂′
 probability at least 3/4. √
√ with
is an O ǫ′′ -additive approximation of p′ (note
√
√
√
√
√
√
that a − γ ≤ a − γ ≤ a + γ ≤ a + γ for
all a ≥ γ ≥ 0). As a result,
p
p
p
(90)
|h1, i|z′ i| − α p̂′ = α p′ − α p̂′
 √ 
= O α ǫ′′
(91)

k=−K

U := i

(87)

2

(76)

So in order to estimate |β̂i | up to additive error O(ǫ),
we only need to obtain an O(ǫ)-additive approximation of |h1, i|z′ i|. This can be achieved as follows.
Let V be a unitary operator such that
J−1 K
1 X X p
|α(j, k)||j, ki,
V |0m i = √
α j=0

,

where W acts on the first and third registers, R acts
on the second and third registers, and Ξ⊥ is an unnormalized state satisfying (|0m ih0m | ⊗ I) Ξ⊥ = 0.
If we measure the first m + 1 qubits of this state in
the standard basis, the probability of getting outcome 0m 1 is

for all j, k.
Now let |zi := A+ |bi = |1i|β̂i and |z′ i :=
h(A)|bi. Then k|zi − |z′ ik = O(ǫ) by Eq. (73).
Thus, for any i ∈ {1, 2, . . . , d}, we have
h1, i|z′ i − β̂i = |h1, i|z′ i − h1, i|zi|

1,2,3

|0m i1 |0i2 |0, ji3
10

Recall that |zi = A+ |bi = |1i|β̂i and |z′ i =
h(A)|bi satisfy √
k|zi − |z′ ik = O(ǫ). As a result,
′
|h1, −i,j |z i| = α p′′ is an O(ǫ)-additive approxima√
tion of |h1, −i,j |zi| = |β̂i − β̂j |/ 2. So in order to
estimate |β̂i − β̂j | up to additive error O(ǫ), we only
need
√ to obtain an O(ǫ)-additive approximation of
α p′′ . To achieve this, we use amplitude estimation to obtain an ǫ′′ -additive
approximation p̂′′ of

′′
′′
2
2
p , where ǫ = Θ ǫ /α , succeeding with proba√
bility at least 3/4.√Then p̂′′ is an O(ǫ/α)-additive
√ ′′
approximation of p′′ , and hence
√ ′′ α p̂ is an O(ǫ)additive approximation of α p , as desired.
The above basic algorithm has success probability
at least 3/4. To raise the success probability to at
least 1 − δ, we repeat this algorithm O(log(1/δ))
times, and take the median of the estimates from
these runs. A standard Chernoff’s bound ensures
that the failure probability is at most δ.
To analyze the complexity of this algorithm, note
that all the parameters are on the same order as in
the proof of Lemma 3. Moreover, Q can be implemented in time poly(log(N )). Therefore, this algorithm makes
 1.5 3
  κ 
d κ
O
(102)
· poly log
ǫ2
ǫδ

implemented in time poly(log(N )). As a result, this
algorithm makes
 1.5 3
  κ 
d κ
O
(94)
· poly log
ǫ2
ǫδ
uses of Px and Py , and is gate-efficient, as claimed.

Lemma 4. Let X, y and β̂ be defined as in LR-P.
Then there exists a gate-efficient quantum algorithm
that makes
 1.5 3
  κ 
d κ
O
· poly log
ǫ2
ǫδ
uses of Px and Py , and outputs an ǫ-additive
approximation of |β̂i − β̂j |, for any given i, j ∈
{1, 2, . . . , d}, succeeding with probability at least 1−δ.
Proof. Let us use the same notation as in the proof of
Lemma 3. The proof of this lemma is quite similar
to that one. The main difference is that here we
replace R with a unitary operator Q satisfying
Q|0i|1, −i,j i = |1i|1, −i,j i,
Q|0i|1, +i,j i = |0i|1, +i,j i,
Q|0i|1, li = |0i|1, li, 1 ≤ l ≤ N, l 6= i, j,
Q|0i|0, ki = |0i|0, ki, 1 ≤ k ≤ N,

(95)
(96)
(97)
(98)

uses of PA and Pb , and is gate-efficient, as claimed.

where
|ii ± |ji
.
|1, ±i,j i := |1i ⊗ √
2

Our algorithm for solving the LR-P problem also
requires the following procedure for determining
whether a given vector β ∈ Rd is close to β̂ or −β̂,
under the promise that one of these cases holds.

(99)

Then by Eqs. (83), (95), (96), (97) and (98), we get

Lemma 5. Let X, y and β̂ be defined as in LR-P.
Suppose β ∈ Rd is given such that either√kβ − β̂k ≤
δ or kβ + β̂k ≤ δ, for some δ < τ /(2σρ d), where
τ := τ (X, y), σ := σ(X), and ρ := ρ(y). Then there
exits a gate-efficient
quantum algorithm that makes

O d1.5 κ uses of Px and Py , and determines which
case holds, succeeding with high probability (e.g. at
least 3/4).

h1, −i,j |z′ i m
|0 i1 |1i2 |1, −i,j i3
α
h1, +i,j |z′ i m
|0 i1 |0i2 |1, +i,j i3
+
α
X h1, l|z′ i
|0m i1 |0i2 |1, li3
+
α

QW |0m i1 |0i2 |bi3 =

l6=i,j

+

X h0, k|z′ i
α

k

⊥

+ Ξ

1,2,3

|0m i1 |0i2 |0, ki3

,

Proof. Let ŷ := Π(X)y = Xβ̂. Then since τ =
2
2
2
τ (X, y) = kŷk /kyk
√ = kŷk (recall that kyk =
1), we have kŷk = τ . Meanwhile, recall that the
singular values of X are in the range [1/κ, 1]. Thus,
we have
√
√
τ ≤ kβ̂k = X+ y = X+ ŷ ≤ κ τ . (103)
√
Note that kβ̂ − (−β̂)k = 2kβ̂k ≥ 2 τ . So by the
√
triangle inequality, at least one of kβ − β̂k ≥ τ
√
and kβ − (−β̂)k ≥ τ must hold. Then since δ <

(100)

where W acts on the first and third registers, Q acts
on the second and third registers, and Ξ⊥ is an
unnormalized state satisfying (|0m ih0m | ⊗ I) Ξ⊥ =
0. If we measure the first m + 1 qubits of this state,
then the probability of getting outcome 0m 1 is
2

p′′ :=

|h1, −i,j |z′ i|
.
α2

(101)
11

√
σρκ τ d
τ
≤
+
N√
2N
2σρκ d
,
≤
N

√
√
τ /(2σρ d) ≤ τ (note that τ ≤ 1 and σ, ρ, d ≥ 1),
the two cases kβ − β̂k ≤ δ and kβ + β̂k ≤ δ cannot
happen simultaneously.
Recall that we have shown in the proof of Lemma
1 that
√
σ d
kxi k ≤ √ , 1 ≤ i ≤ N
(104)
N

xTi β̂

√
σκ τ d
, 1 ≤ i ≤ N.
≤ √
N

N
X
i=1

(105)

=

i=1

N
X
i=1
T

yi · xTi β̂

= y Xβ̂
= yT Π(X)y
= kŷk
= τ.

2

≤ |yi |kxi kkβ − β̂k
√
σ d
ρ
≤ √ · √ ·δ
N
N
τ
<
,
2N

(108)
(109)
(110)

τ
.
2

τ
2

i=1

|qi − q̂i |

(120)
(121)
(122)

(124)

(125)
(126)

for all i ∈ {1, 2, . . . , N }. Then by Eqs. (107)
and (126), we find that

(111)
(112)

|qi | < |q̂i | + |qi + q̂i |
√
σρκ τ d
τ
≤
+
N√
2N
2σρκ d
≤
,
N

(127)
(128)
(129)

for all i ∈ {1, 2, . . . , N } (note that ρ, σ, κ, d ≥ 1
and τ ≤ 1). Furthermore, Eqs. (112) and (126)
imply that
N
X

(113)
(114)

i=1

qi ≤ −

N
X

τ
=− .
2

(116)

q̂i +

N
X
i=1

i=1

< −τ +

(115)

τ
2

|qi + q̂i |

(130)
(131)
(132)

Comparing Eqs. (122) and (132), we know that we
can distinguish the two cases kβ − β̂k ≤ δ and
P
kβ + β̂k ≤ δ by estimating N
i=1 qi up to additive
error τ /2, as claimed.

for all i ∈ {1, 2, . . . , N }. Then by Eqs. (107)
and (116), we find that
|qi | < |q̂i | + |qi − q̂i |

q̂i −

≤ |yi |kxi kkβ + β̂k
√
σ d
ρ
≤ √ · √ ·δ
N
N
τ
,
<
2N

Now let qi := yi · xTi β for i ∈ {1, 2, . . . , N }. We
claim that we can distinguish the cases kβ − β̂k ≤ δ
PN
and kβ + β̂k ≤ δ by estimating the quantity i=1 qi
up to additive error τ /2. To prove this, let us consider these two cases separately:
√
• Case 1: kβ − β̂k ≤ δ < τ /(2σρ d). Using
Eqs. (104) and (106), we get


|qi − q̂i | = yi · xTi β − β̂

i=1

N
X

√
• Case 2: kβ + β̂k ≤ δ < τ /(2σρ d). Using
Eqs. (104) and (106), we get


|qi + q̂i | = yi · xTi β + β̂
(123)

Now let q̂i := yi · xTi β̂ for i ∈ {1, 2, . . . , N }. Then
Eqs. (105) and (106) imply that
√
σρκ τ d
, 1 ≤ i ≤ N,
(107)
|q̂i | ≤
N
Furthermore, we have
q̂i =

qi ≥

N
X

>τ−

Moreover, by kyk = 1 and ρ(y) = ρ, we obtain
ρ
|yi | ≤ √ , 1 ≤ i ≤ N.
(106)
N

N
X

(119)

for all i ∈ {1, 2, . . . , N } (note that ρ, σ, κ, d ≥ 1
and τ ≤ 1). Furthermore, Eqs. (112) and (116)
imply that

Combining Eqs. (103) and (104)

(see Eq. (47)).
yields

(118)

(117)
12

Theorem 2. The LR-P problem can be solved by a
gate-efficient quantum algorithm that makes
 2.5 3
  
d κ
dκ
O
· poly log
2
δ
δ

We obtain a τ /2-additive approximation of
PN
i=1 qi as follows. Let U be a unitary operator such
that
U |ii|0i = |ii|ψi i, 1 ≤ i ≤ N,

(133)

where

uses of Px and Py , where δ := min{ǫ, 1/d}.

r

r
1 N qi
1 N qi
|ψi i :=
+
|0i +
−
|1i (134)
2
2∆
2
2∆
√
in which ∆ := 2σρκ d. Note that U is a valid
unitary operator, since N |qi | ≤ ∆ by Eqs. (119) and
(129). Then we have
!
N
N
1 X
1 X
U √
|ii |0i = √
|ii|ψi i. (135)
N i=1
N i=1

Proof. Algorithm: Let X and y be defined as in
LR-P. Let τ := τ (X, y) = Ω(1), σ := σ(X) = O(1)
and ρ := ρ(y) = O(1). We use the following algorithm to obtain a vector β = (β1 , β2 , . . . , βd )T ∈ Rd
satisfying kβ − β̂k∞ ≤ ǫ, succeeding with probability
at least 2/3:
1. Let ǫ′ := min{τ /(2σρd), ǫ}.
2. For each j ∈ {1, 2, . . . , d}, we run the algorithm in Lemma 3 to obtain an ǫ′ /6-additive
approximation µj of |β̂j |, succeeding with
probability at least 1 − 1/(25d).

If we measure the second register of this state in
the standard basis, then the probability of obtaining
outcome 0 is
PN
qi
1
(136)
p := + i=1 .
2
∆

3. Let S := {j ∈ {1, 2, . . . , d} : µj > 2ǫ′ /3}. If
S = Ø, then this algorithm fails; otherwise,
we continue as follows.

We use amplitude estimation to obtain an τ /(2∆)additive approximation p̂ of p, succeeding with high
probability (e.g. at least 3/4). Then (p̂ − 1/2)∆ is a
PN
τ /2-additive approximation of i=1 qi , as desired.
The unitary operator U can be implemented as
follows. For any i ∈ {1, 2, . . . , N }, given the state
|ii|0i, we first transform it into |ii|0i|qi i, where qi =
Pd
yi ( j=1 xi,j βj ) can be computed by making O(d)
uses of Px and Py . Then we perform the controlledrotation
|0i|qi i → |ψi i|qi i

4. Pick arbitrary j0 ∈ S. For each j ∈ S,
j 6= j0 , we run the algorithm in Lemma 4
to obtain an ǫ′ /6-additive approximation γj of
|β̂j0 − β̂j |, succeeding with probability at least
1 − 1/(25d).
5. For each j ∈ {1, 2, . . . , d}, we define sj ∈
{−1, 0, 1} as follows:
• If j 6∈ S, then sj = 0.

(137)

• If j = j0 ∈ S, then sj = 1.

on the last two registers. After that, we uncompute
qi in the last register by making O(d) uses of Px and
Py , and get the desired state |ii|ψi i. This implemenation of U requires O(d) uses of Px and Py , and is
gate-efficient.
Since we want to estimate p up to additive error
τ /(2∆), amplitude estimation requires
√ !
 
 √ 
σρκ d
∆
=O
=O κ d
O
(138)
τ
τ

• Otherwise, we have j ∈ S and j 6= j0 .
If ||µj0 − µj | − γj | ≤ ǫ′ /2, then sj = 1;
otherwise, sj = −1.

repetitions of U (recall that σ = O(1), ρ = O(1)
and τ = Ω(1)). As a result, this algorithm makes
O d1.5 κ uses of Px and Py , and is gate-efficient, as
claimed.

6. Let β ′ = (β1′ , β2′ , . . . , βd′ )T ∈ Rd be defined
as βj′ := sj µj for each j ∈ {1, 2, . . . , d}.
We will prove below that, with high√ probability, either kβ ′√− β̂k < τ /(2σρ d) or
kβ ′ + β̂k < τ /(2σρ d). We run the algorithm
in Lemma 5 to determine which case holds,
succeeding with probability at least 3/4. If
the first case holds, then we return β := β ′
as our estimate of β̂; otherwise, we return
β := −β ′ as our estimate of β̂.

Now we are ready to state our algorithm for solving the LR-P problem.

Correctness: Let us call the case where all
the instances of the algorithms in Lemmas 3, 4
13

and 5 succeed the typical case. By union bound,
the probability of this case happening is at least
1 − 2d/(25d) − 1/4 > 2/3. We will prove that in
the typical case, our algorithm outputs a correct β
(i.e. kβ − β̂k∞ ≤ ǫ) with certainty.
In the typical case, we have
µi − β̂i ≤

ǫ′
, 1 ≤ i ≤ d,
6

• If j = j0 , then sj = 1 by definition.
• If j 6= j0 and sgn(β̂j ) = sgn(β̂j0 ), then we have
β̂j0 − β̂j = β̂j0 − β̂j .

Combining Eqs. (139), (140) and (152) gives

(139)

||µj0 − µj | − γj | ≤ µj0 − β̂j0

and
γj − β̂j0 − β̂j

≤

ǫ′
, ∀j ∈ S, j 6= j0 .(140)
6

+ µj − β̂j
+ γj − β̂j0 − β̂j

Then using the definition of S, we get
β̂j ≥ |µj | − µi − β̂i
ǫ′
2ǫ′
−
3
6
ǫ′
= , ∀j ∈ S,
2
>

(152)

′

′

(153)

′

ǫ
ǫ
ǫ
+ +
6
6
6
ǫ′
= .
2

≤

(141)
(142)

(154)
(155)

This implies that sj = 1 for this j.

(143)

• If j 6= j0 and sgn(β̂j ) = −sgn(β̂j0 ), then we
have

and
β̂j ≤ |µj | + µi − β̂i
2ǫ′
ǫ′
≤
+
3
6
5ǫ′
, ∀j 6∈ S.
=
6

(144)
β̂j0 − β̂j = β̂j0 + β̂j

(145)

= β̂j0 − β̂j
n
o
+ 2min β̂j0 , β̂j

(146)

Recall that we have shown in the proof of Lemma
5 that
v
u d
uX
√
β̂ = t
(147)
|β̂ |2 ≥ τ

(156)

> β̂j0 − β̂j + ǫ′ ,

(157)

since |β̂j0 |, |β̂j | > ǫ′ /2 by Eq. (143). Combining Eqs. (139), (140) and (157) yields

i

i=1

||µj0 − µj | − γj | ≥ ǫ′ − µj0 − β̂j0

(see Eq. (103)). This implies that there exists some
i0 ∈ {1, 2, . . . , d} such that
r
τ
τ
β̂i0 ≥
≥
≥ 2ǫ′
(148)
d
σρd

− µj − β̂j
− γj − β̂j0 − β̂j

(note that σ, ρ, d ≥ 1 and τ ≤ 1). Then by Eqs. (139)
and (148), we obtain
µi0 ≥ |βi0 | − µi0 − β̂i0
≥ 2ǫ′ −
>

2ǫ′
.
3

ǫ′
6

> ǫ′ −
=

(149)

ǫ′
.
2

′

′

(158)
′

ǫ
ǫ
ǫ
− −
6
6
6

(159)
(160)

This implies that sj = −1 for this j.

(150)

The fact that sj = sgn(β̂j ) · sgn(β̂j0 ) for all j ∈ S
implies that either

(151)

sgn(βj′ ) = sgn(β̂j ),

Thus, we have i0 ∈ S and S 6= Ø. So our algorithm
does not fail in the typical case.
Now we claim that sj = sgn(β̂j ) · sgn(β̂j0 ) for any
j ∈ S. The proof is as follows.

∀j ∈ S,

(161)

or
sgn(βj′ ) = −sgn(β̂j ),
14

∀j ∈ S.

(162)

Moreover, by Eq. (139), we know that
ǫ′
|βj′ | − |β̂j | ≤ ,
6
As a result, we have either
βj′ − β̂j ≤

ǫ′
,
6

∀j ∈ S.

∀j ∈ S,

=O
(163)

(164)

ǫ′
, ∀j ∈ S.
(165)
6
Meanwhile, for any j 6∈ S, we have sj = 0 and
|β̂j | ≤ 5ǫ′ /6 by Eq. (146). It follows that βj′ = 0 and
βj′ + β̂j ≤

Our algorithm for computing β̂ = X+ y is more
efficient than an alternative one in which one first
creates multiple copies of the state proportional to
β̂ and then uses statistical sampling and quantum
state tomography to determine the β̂j ’s (as suggested by Ref. [1]). The main reason is that, in order to obtain an ǫ-additive approximation of |β̂j |2,
the sampling-based approach would require O 1/ǫ2
copies of the state encoding β̂, but amplitude estimation only needs O(1/ǫ) repetitions of the procedure
for preparing this state. So it is more efficient to
couple the state generation process with amplitude
estimation (as we did in our algorithm) rather than
statistical sampling.
We also remark that the algorithm in Lemma 3
can be modified to produce a quantum state approximately proportional to β̂. Specifically, note
that if we measure the first register of W |0m i|bi (in
Eq. (83)) in the standard basis, then conditioning on
the outcome being 0m , we would obtain the normalized version of h(A)|bi, which is close to the normalized version of A+ |bi = |1i|β̂i. The probablity of

2
this event happening is kh(A)|bik /α2 = Ω 1/α2 .
We can use amplitude amplification to raise this
probability to Ω(1), which requires O(α) repetitions
of W and the procedure for preparing |bi. This leads
to a gate-efficient algorithm that makes

  κ 
(176)
O d1.5 κ2 · poly log
ǫ
uses of Px and Py , and prepares a quantum state

5ǫ′
, ∀j 6∈ S. (166)
+ β̂j ≤
− β̂j =
6
Combining the cases j ∈ S and j 6∈ S, we know
that either
5ǫ′
< ǫ′
6

(167)

kβ ′ + β̂k∞ ≤

5ǫ′
< ǫ′ .
6

(168)

or

As a result, we have either
√
kβ ′ − β̂k < dǫ′ ≤
or
kβ ′ + β̂k <

√ ′
dǫ ≤

τ
√
2σρ d

(169)

τ
√ .
2σρ d

(170)

(174)

uses of Px and Py , and is gate-efficient, as claimed.

βj′

kβ ′ − β̂k∞ ≤

  
d2.5 κ3
dκ
·
poly
log
δ2
δ

uses of Px and Py , and is gate-efficient. Step 6 makes
O(1) uses
 of the algorithm in Lemma 5, so it requires
O κd1.5 uses of of Px and Py , and is gate-efficient.
Furthermore, the classical computation in this algorithm takes O(d) time. As a result, this algorithm
makes
 2.5 3
  
d κ
dκ
O
(175)
· poly log
2
δ
δ

or

βj′



In the typical case, our algorithm in Lemma 5
correctly determines which case holds. If the first
case holds, then it outputs β = β ′ which satisfies
kβ − β̂k∞ < ǫ′ ≤ ǫ; otherwise, it outputs β = −β ′
which also satisfies kβ − β̂k∞ < ǫ′ ≤ ǫ, as desired.
Complexity: Recall that ǫ′ = min{τ /(2σρd), ǫ}
and δ = min{1/d, ǫ}, where τ = Ω(1), σ = O(1) and
ρ = O(1). So we have ǫ′ = Ω(δ). Let us analyze the
complexity of each step. Step 2 makes O(d) uses of
the algorithm in Lemma 3, so it requires

  
d1.5 κ3
dκ
O d·
(171)
·
poly
log ′
(ǫ′ )2
ǫ
 2.5 3
  
d κ
dκ
=O
(172)
· poly log
δ2
δ

|β̂i
in l2 norm, succeeding with probak|β̂ik
bility Ω(1) (with a flag indicating success). By utilizing Ambainis’ variable-time amplitude amplification
[24], we can reduce the κ-dependence from quadratic
to linear, as done in Section 5 of Ref. [7]. This leads
to a gate-efficient algorithm with query complexity

  κ 
(177)
O d1.5 κ · poly log
ǫ

ǫ-close to

uses of Px and Py , and is gate-efficient. Step 4 makes
O(d) uses of the algorithm in Lemma 4, so it requires
  

dκ
d1.5 κ3
(173)
· poly log ′
O d·
(ǫ′ )2
ǫ
15

Proof. We can get a ∆/2-additive approximation of
θj by using the standard phase estimation, which
makes O(1/∆) uses of U and uses poly(log(1/∆))
additional 2-qubit gates. This is sufficient to distinguish the two cases. However, it only succeeds with
probability Ω(1). To raise this probability to at least
1−δ, we repeat this procedure O(log(1/δ)) times and
check whether the median of the estimates is ∆/2close to θ. A standard Chernoff’s bound ensures
that the failure proability is at most δ. This boosted
procedure, denoted by P, makes O((1/∆) · log(1/δ))
uses of U , and uses poly(log(1/(∆δ))) additional
2-qubit gates, and satisfies all the desired properties.

for the same task.
One may compare this algorithm for preparing
a quantum state approximately proportional to the
optimal parameters
β̂ = X+ y = (XT X)−1 XT y

(178)

with the one in Ref. [1] for the same task. Our algorithm is based on the singular value decomposition
(SVD) of X, and it applies X+ to y in a direct manner. Consequently, it has only linear dependence on
the condition number κ of X. By contrast, Ref. [1]
needs to first apply XT to y, which incurs a κ factor
in the complexity; then it needs to apply (XT X)−1
to the output of the first step, which incurs another
κ2 factor in the complexity. So its overall complexity
is at least cubic in κ. This means that our algorithm
has polynomially better dependence on κ than the
one in Ref. [1]. Furthermore, due to the fact we use
the new strategy of Ref. [7] for matrix inversion, our
algorithm also has exponential better dependence on
the desired precision ǫ in the output state.
V.

Theorem 3. The LR-Q problem can be solved by a
gate-efficient quantum algorithm that makes
 1.5
  κ 
d κ
O
· poly log
ǫ
ǫ
uses of Px and Py .
Proof. Algorithm: Let X and y be defined as in
LR-Q. We use the following algorithm to obtain an
ǫ-additive approximation of τ = kΠ(X)yk2 /kyk2 =
2
kΠ(X)yk (recall that kyk = 1), succeeding with
probability at least 2/3. Let A := |1ih0| ⊗ XT +
|0ih1| ⊗ X and |bi := |0i|yi. Let P be the unitary
procedure in Lemma 6 for U = e−iA , θ = 0, ∆ =
1/(2κ) and δ = ǫ/2. Suppose

ESTIMATING THE QUALITY OF THE
LEAST-SQUARES FIT

In this section, we describe a quantum algorithm
for solving the LR-Q problem, i.e. estimating the
quality τ = kXβ̂k2 /kyk2 of the least-squares fit
y ≈ Xβ̂ for a given data set (X, y) (without computing the parameters β̂ explicitly). This algorithm
requires the following variant of phase estimation
[25, 26], which decides whether the eigenphase corresponding to an eigenvector of a unitary operator
is θ or far away from θ, for some given θ ∈ [0, 2π),
succeeding with probability close to 1. (Similar procedures have been used in Refs. [7, 27, 28].)

P|0i1 0l 2 |bi3 = µ0 |0i1 |ϕ0 i2,3 + µ1 |1i1 |ϕ1 i2,3 ,
(180)
2
2
where l = O(log(1/∆) log(1/δ)), |µ0 | + |µ1 | = 1,
and |ϕ0 i2,3 and |ϕ1 i2,3 are some normalized states
on the second and third registers. We use amplitude
estimation to get an ǫ/2-additive approximation r̂
2
of r := |µ1 | , succeeding with probability at least
3/4. Then we return r̂ as our estimate of τ . During
this process, we use the procedure in Lemma 1
to implement U = e−iA with precision
O ǫ2 /κ2

2
2
(and failure probability O ǫ /κ ), and use the
procedure
in Lemma 2 to prepare |yi with precision

O ǫ2 .

Lemma 6. Let U be a unitary operator with eigenvectors |ψj i such that U |ψj i = eiθj |ψj i for some
θj ∈ [0, 2π). Let θ ∈ [0, 2π) and let ∆, δ ∈
(0, 1). Then there is a unitary procedure P that
makes O((1/∆) · log(1/δ)) uses of U , and uses
poly(log(1/(∆δ))) additional 2-qubit gates, and satisfies

Correctness: Suppose X has the singular value decomposition

P|0i 0l |ψj i = (αj,0 |0i|ηj,0 i + αj,1 |1i|ηj,1 i)|ψ(179)
j i,
2

2

where l = O(log(1/∆) log(1/δ)), |αj,0 | +|αj,1 | = 1,
|ηj,0 i and |ηj,1 i are two normalized states, and

X=

d
X
j=1

2

• If θj = θ, then |αj,0 | ≥ 1 − δ.

sj |uj ihvj |,

(181)

where sj ∈ [1/κ, 1], |uj i ∈ RN and |vj i ∈ Rd are
unit vectors, for all j ∈ {1, 2, . . . , d}. Then A has

2

• If |θj − θ| ≥ ∆, then |αj,1 | ≥ 1 − δ.
16

the spectral decomposition
A=

d
X
j=1

sj |+j ih+j | −

d
X
j=1

+
sj |−j ih−j |, (182)

d
X

αj −
−
−
√ γj,0
|0i φ−
j,0 + γj,1 |1i φj,1 |−j i
2
j=1

+ α(η0 |0i|ψ0 i + η1 |1i|ψ1 i)|0i Φ⊥ .

where

(192)

It follows that
1
|±j i := √ (|0i|uj i ± |1i|vj i).
2

(183)

d

r=

Meanwhile, we can write |yi as
|yi =
Pd

2

d
X
j=1

αj |uj i + α Φ⊥ ,
2

2

d
X
j=1

2

|αj | .

|r − τ | ≤

(185)
+

By Eqs. (183) and (184), we obtain
|bi = |0i|yi
=

d
X
j=1

=

≤

(187)

j=1



2

2

+ |α| |η1 | . (193)

j=1

2

(194)
2

|αj | · δ +

+ |α| · δ
=δ
ǫ
= .
2

1
2

d
X
j=1

2

|αj | · δ
(195)
(196)
(197)

Namely, r is an ǫ/2-additive approximation of τ .
Meanwhile, r̂ is an ǫ/2-additive approximation of r.
It follows that r̂ is an ǫ-additive approximation of τ ,
as desired.
In the above argument, we have ignored the error
in the implementation of U = e−iA and the error
in the preparation of |yi. We will show below that
our algorithm
only makes o κ2 /ǫ2 uses of U and

2
o 1/ǫ uses of the procedure for preparing |yi.
Thus, provided
 that U is implemented with preci
sion O ǫ2 /κ2 (and failure probability O ǫ2/κ2 )
and |yi is prepared with precision O 1/ǫ2 , the
error in the final state (compared to the ideal case)
is only o(1). Consequently, our algorithm outputs a
correct r̂ (i.e. |r̂ − τ | ≤ ǫ) with probability at least
3/4 − o(1).

P|0i 0l |0i Φ⊥ = (η0 |0i|ψ0 i + η1 |1i|ψ1 i)|0i Φ⊥ ,
(191)
2
2
where |η0 | ≥ 1−δ, |η1 | ≤ δ, |ψ0 i and |ψ1 i are some
normalized states. As a result, we have
αj +
+
+
√ γj,0
|0i φ+
j,0 + γj,1 |1i φj,1
2

1
2

d
X
2

Note that |0i Φ⊥ is an eigenvector of A with eigenvalue 0, i.e. A|0i Φ⊥ = 0.
Now, since the eigenphase gap around 0 of U =
e−iA is at least 1/κ, by Lemma 6 and our choice of
parameters, we get

+
+
+
|0i φ+
P|0i 0l |+j i = γj,0
j,0 + γj,1 |1i φj,1 |+j i,
(189)

−
−
−
|0i φ−
P|0i 0l |−j i = γj,0
j,0 + γj,1 |1i φj,1 |−j i,
(190)
±
± 2
± 2
where γj,1 ≥ 1 − δ, γj,0 ≤ δ, φj,0 and φ±
j,1
are some normalized states, for all j ∈ {1, 2, . . . , d},
and

P|0i 0l |bi =

2

d


1X
2
− 2
|αj | 1 − γj,1
2 j=1
2

d
X
αj
√ (|+j i + |−j i) + α|0i Φ⊥ . (188)
2
j=1

d
X

−
+ γj,1

d


1X
+ 2
|αj |2 1 − γj,1
2 j=1

+ |α| |η1 |

(186)

αj |0i|uj i + α|0i Φ⊥

2

± 2
Note that since |γj,1
| ≈ 1 and |η1 | ≈ 0, we have
r ≈ τ by Eqs. (185) and (193). More precisely, the
difference between r and τ can be bounded using the
triangle inequality:

(184)

where j=1 |αj | + |α| = 1, and Φ⊥ is some normalized state satisfying uj |Φ⊥ = 0 for all j. Note
that
τ = kΠ(X)yk =


1X
2
+
|αj | γj,1
2 j=1

Complexity: Since we want to estimate r up to
additive error O(ǫ), amplitude estimation requires

|+j i O(1/ǫ) repetitions of the procedure P and the procedure for preparing |yi. Then by Lemma 6, our
17

whether z is in the range [0, ⌊N/2⌋] or [⌊N/2⌋+1, N ].
If the first case holds, then we test whether z is in
the range [0, ⌊N/4⌋] or [⌊N/4⌋+1, ⌊N/2⌋]; otherwise,
we test whether z is in the range [⌊N/2⌋+ 1, ⌊3N/4⌋]
or [⌊3N/4⌋ + 1, N ], and so on. We only need
O(log(N )) such tests to locate z, since each test
reduces the size of candidate set by a factor of 2.
Furthermore, by assumption, each test can be accomplished by making at most Q queries to Pf .
Thus, we can find z by making O(Q log(N )) queries
to Pf . On the other hand, it is known that unstruc√ 
tured search has quantum query complexity Ω N

algorithm makes


 
 
1 1
κ
1
1
O
=O
(198)
· log
· log
ǫ ∆
δ
ǫ
ǫ
uses of U . By Lemma 1, U = e−iA can be implemented with precision
O ǫ2 /κ2 (and failure prob
ability O ǫ2 /κ2 ) by a gate-efficient procedure that
makes O d1.5 · log(κ/ǫ) uses of Px . Meanwhile, by
Lemma 2, |yi can be prepared with precision O ǫ2
by a gate-efficient procedure that makes O(log(1/ǫ))
uses of Py . As a result, this algorithm makes
 1.5
  κ 
d κ
(199)
O
· poly log
ǫ
ǫ

[29, 30].
√Combiningthese two facts, we know that
Q = Ω N / log(N ) .

uses of Px and Py , and is gate-efficient, as claimed.

Theorem 4. The LR-P problem has quantum query
complexity Ω(κ/ log(κ)), where κ is the condition
number of the design matrix X.

Comparing Theorem 2 and Theorem 3, one can
see that it is easier to estimate the quality of the
least-squares fit y ≈ Xβ̂ than to find its parameters
β̂ = X+ y explicitly. Thus, in practice, we can first
run the algorithm in Theorem 3 to check whether
a given data set is well-behaved (e.g. τ ≥ 2/3).
If so, then we run the algorithm in Theorem 2 to
fit a linear regression model to this data set. The
total cost of this process is dominated by that of the
second stage.

Proof. We prove this theorem by showing that for
any positive integer N , there exists a balanced
N ×2
matrix
with singular values s1 (X) =
 √ X ∈ R
Θ 1/ N and s2 (X) = Θ(1) such that, for y =
√1 (1, 1, . . . , 1)T
N
T

∈ RN , β̂ = X+ y is either (1, 0)T or

√
(0, 1) , but one has to make Ω N /log(N ) queries
to X to determine which case holds (succeeding with
probability at least 2/3).
Let X √
be an N × 2 matrix such that its entries
are all 1/ N except one entry 0 (whose location is
unknown and arbitrary). Then we know that one
column of X is equal to y = √1N (1, 1, . . . , 1)T , and
the other column of X is linearly independent from
y. Consequently, using the definition

VI. LOWER BOUND ON THE
COMPLEXITY OF LINEAR REGRESSION

Our quantum algorithm for computing β̂ = X+ y
has polynomial dependence on the condition number κ of the design matrix X. In this section, we
show that this dependence is indeed necessary. To
prove this, we need the following lower bound on
the quantum query complexity of a weaker version
of unstructured search.

β̂ = argminkXβ − yk,

(200)

β∈R2

Lemma 7. Let f : {1, 2, . . . , N } → {0, 1} be a function such that f (x) = 1 if and only if x = z for some
unknown z ∈ {1, 2, . . . , N }. Let Pf be a procedure
that on input x ∈ {1, 2, . . . , N }, outputs
√ the value

of f (x). Then one has to make Ω N / log(N )
queries to Pf to determine whether the unknown z
is larger than ⌊N/2⌋ or not (succeeding with probability at least 2/3).

we obtain that β̂ is either (0, 1)T or (1, 0)T , depending on whether the entry 0 is in the first or second
column 
of X, respectively.
By Lemma 7, one must

√
make Ω N / log(N ) queries to X to determine
which column contains theentry 0. This implies that
√
one also needs to make Ω N / log(N ) queries to
X to determine whether β̂ = (0, 1)T or β̂ = (1, 0)T .
It remains to show that X also satisfies the other
desired properties. First, by a direct calculation,
 √  we
get that kXkF = Θ(1), kXk2,∞ = Θ 1/ N and

Proof. Suppose we can solve the given problem by
making Q queries to Pf . Then we can find the unknown z by making O(Q log(N )) queries to Pf . The
idea is to use binary search. Namely, we first test
18

produce a quantum state encoding the optimal parameters, our algorithm outputs these numbers in
the classical form. So by running it once, one completely determines the fitted model and then can use
it to make predictions on new data at little cost.
The running time of this algorithm is polynomial in
log(N ), d, κ and 1/ǫ, where N is the size of the data
set, d is the number of adjustable parameters, κ is
the condition number of the design matrix, and ǫ is
the desired precision in the output. We also show
that the polynomial dependence on d and κ is necessary. Therefore, our algorithm cannot be greatly
improved. Furthermore, we also give an efficient
quantum algorithm that estimates the quality of the
least-squares fit (without computing its parameters
explicitly). This algorithm runs faster than the one
for finding this fit, and can be used to check whether
the given data set qualifies for linear regression in the
first place.
One may have noticed that our algorithms actually solve two fundamental problems in linear algebra. One is to apply the pseudoinverse of a dense
rectangular matrix to a vector, and the other is to
estimate the norm of the projection of this vector
onto the range of this matrix. Such problems frequently arise in many scenarios. So it is conceivable
that our algorithms may find applications beyond
linear regression.
Our algorithms might be improved in a few ways.
Ambainis [24] proposed a technique called variabletime amplitude amplification and utilized it to enhance the κ-dependence of HHL’s algorithm [2] for
preparing a state encoding the solution of a linear
system (this techique is also used in CKS’s algorithm
[7]). But it is unknown whether this technique leads
to a more efficient algorithm for estimating an entry
(or the difference between two entries) of this solution. If so, we would obtain a faster algorithm for
fitting a linear regression model to a data set using
the least squares approach. On the other hand, for
estimating the quality of the fitted model, we still
do not know whether the polynomial dependence on
κ is necessary. We believe that this is the case, but
could not prove it. This is left as an interesting open
question.
In this paper, we have focused on linear regression with ordinary least squares optimization (which
assumes that the errors for different observations
are independent). It is also worth investigating the
quantum complexity of linear regression with generalized least squares optimization (which allows the
errors for different observations to be correlated).
Furthermore, one might study how these complexities change when regularization is used. For ex-

hence σ(X) = Θ(1). Second, note that either


1
1
1 − N 1 − N 
(201)
XT X = 

1
1
1−
N

or


1
1− 
 1
N .
XT X = 
1
1
1−
1−
N
N


(202)

By a direct calculation, we find that λ1 (XT X) =
Θ(1/N ) and λ2 (XT X) = Θ(1). It follows that
p
√
s1 (X) =
λ1 (XT X) = Θ 1/ N and s2 (X) =
√ 
p
λ2 (XT X) = Θ(1), and hence κ(X) = Θ N .
This concludes the proof.
Clearly, the LR-P problem has time complexity
Ω(d), because simply writing down a d-dimensional
vector β ≈ β̂ requires this amount of time. Combining this fact and Theorem 4, we know that the
algorithm in Theorem 2 cannot be dramatically improved.
It is worth noting that Harrow, Hassidim and
Lloyd (HHL) [2] have also given a lower bound on
the quantum complexity of matrix inversion. They
proved that unless BQP = PSPACE, one cannot
solve the matrix inversion problem in quantum time
κ1−δ · poly(log(N )) for some constant δ > 0, where
κ and N are the condition number and dimension of
the matrix to be inverted, respectively. We remark
that this result and Theorem 4 are incomparable. At
first glance, it may seem that Theorem 4 is stronger,
since it has better dependence on κ and it does not
rely on any complexity-theoretic assumption. But
recall that in our LR-P problem, we allow the design
matrix X to be nonspare, while HHL only allowed
sparse matrices in their work. So we only obtain a
stronger bound under a stronger assumption. Nevertheless, it may be possible to use our approach
to improve HHL’s bound, showing that our bound
holds in the sparse case as well. This is left as an
interesting open question.

VII.

DISCUSSION

To summarize, we have presented an efficient
quantum algorithm for fitting a linear regression
model to a given data set using the least squares
approach. Different from previous algorithms which
19

ample, how hard is it to solve ridge regression [31]
or Lasso [32] on a quantum computer? Finally, it
would be worth exploring the power and limitation
of quantum algorithms for nonlinear regression.
Our work is also a new contribution to the nascent
field of quantum machine learning, which has made a
lot of progress in the past years [1, 3, 33–57]. Here we
briefly review this broad area, and position our work
with the other works in this area (for an excellent
review on quantum machine learning, see Ref. [50]).
In fact, depending on the types of the learning device and the object to be learned, quantum machine
learing can be divided into three branches. The first
branch, which is also known as quantum-enhanced
machine learning, uses quantum mechanics to improve the performance of classical machine learning
methods (e.g. [1, 3, 41, 42, 52]). Conversely, the second branch applies classial machine learning methods to the study of quantum systems (e.g. [33, 43]).
Finally, the third branch uses quantum approaches
to study quantum systems (e.g. [46]). Clearly, our
work is an instance of the first kind, i.e. quantumenhanced machine learning.
Now let us look at quantum-enhanced machine
learning more carefully. Traditional machine learning algorithms can be divided into three main groups
based on their purpose: supervised learning (in
which an algorithm learns from example data and
associated target responses that can consist of numeric values or string labels), unsupervised learning (in which an algorithm learns from plain examples without any associated response), and reinforcement learning (in which an agent interacts with an
environment and occasionally receives rewards for
its actions, which allows the agent to adapt its behavior). There has beening exciting progress in all
of these three paradigms. See Refs. [1, 3, 42, 47],
Refs. [36, 37, 52] and Refs. [34, 41, 53, 54] for examples of the first, second and third kind, respectively.
Our work belongs to the first category, as it concerns
least-square linear regression – a typical supervised
learning task.

Meanwhile, we can also classify the works on
quantum-enhanced machine learning based on the
techniques they use. It seems that most of these
works fall into three groups according to this criterion. The first group use linear algebra methods (e.g.
singular value decomposition), and are usually related to HHL’s quantum algorithm for linear systems
of equations somehow. Examples include Refs. [1, 3]
and this work on least-squares linear regression,
Ref. [42] on support vector machine, and Ref. [47]
on Guassian processes. This approach could achieve
exponential speedup (in some sense) over classical
methods. The second group are based on amplitude amplification (including Grover’s search and
quantum walk). Examples include Ref. [37] on kmedians, Ref. [49] on k-nearest neighbors, Ref. [35]
on Google’s PageRank, and Ref. [41] on reinforcement learning. This approach usually achieves polynomial speedup over classial methods. Finally, the
third group are based on quantum sampling techniques (e.g. quantum annealing). Examples include
Refs. [48, 52, 55–57] on (deep) Boltzmann machines.
We believe that the field of quantum(-enhanced) machine learning could benefit the most from the marriage of these different ideas, and look forward to
seeing more novel quantum algorithms for solving
machine learning tasks.

[1] N. Wiebe, D. Braun, and S. Lloyd, Phys. Rev. Lett.
109, 050505 (2012).
[2] A. W. Harrow, A. Hassidim, and S. Lloyd, Phys.
Rev. Lett. 103, 150502 (2009).
[3] M. Schuld, I. Sinayskiy, and F. Petruccione, Phys.
Rev. A 94, 022342 (2016).
[4] S. Lloyd, M. Mohseni, and P. Rebentrost, Nature
Physics 10, 631 (2014).

[5] G. H. Low and I. L. Chuang, arXiv preprint
arXiv:1610.06546 (2016).
[6] G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118,
010501 (2017).
[7] A. M. Childs, R. Kothari, and R. D. Somma, arXiv
preprint arXiv:1511.02306 (2015).
[8] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp,
Contemporary Mathematics 305, 53 (2002).

ACKNOWLEDGMENTS

The author thanks Scott Aaronson, Andrew
Childs and Umesh Vazirani for helpful discussions
and comments. The author also thanks the anonymous referee for providing many useful comments
on an earlier version of this paper. Part of this work
was done while the author was a graduate student
at Computer Science Division, University of California, Berkeley. This research was supported by ARO
Grant W911NF- 09-1-0440.

20

[9] G. H. Golub and C. Reinsch, Numerische mathematik 14, 403 (1970).
[10] S. Lloyd, Science 273, 1073 (1996).
[11] D. Aharonov and A. Ta-Shma, in Proceedings of the
Thirty-fifth Annual ACM Symposium on Theory of
Computing, STOC ’03 (ACM, New York, NY, USA,
2003) pp. 20–29.
[12] A. M. Childs, Quantum information processing in
continuous time, Ph.D. thesis, Massachusetts Institute of Technology (2004).
[13] D. W. Berry, G. Ahokas, R. Cleve, and B. C.
Sanders, Communications in Mathematical Physics
270, 359 (2007).
[14] A. M. Childs, Communications in Mathematical
Physics 294, 581 (2010).
[15] A. M. Childs and R. Kothari, “Simulating sparse
hamiltonians with star decompositions,” in Theory of Quantum Computation, Communication, and
Cryptography: 5th Conference, TQC 2010, Leeds,
UK, April 13-15, 2010, Revised Selected Papers,
edited by W. van Dam, V. M. Kendon, and S. Severini (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 94–103.
[16] D. Poulin, A. Qarry, R. Somma, and F. Verstraete,
Phys. Rev. Lett. 106, 170501 (2011).
[17] A. M. Childs and N. Wiebe, Quantum Info. Comput. 12, 901 (2012).
[18] D. W. Berry and A. M. Childs, Quantum Info. Comput. 12, 29 (2012).
[19] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari,
and R. D. Somma, in Proceedings of the Forty-sixth
Annual ACM Symposium on Theory of Computing,
STOC ’14 (ACM, New York, NY, USA, 2014) pp.
283–292.
[20] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari,
and R. D. Somma, Phys. Rev. Lett. 114, 090502
(2015).
[21] D. W. Berry, A. M. Childs, and R. Kothari, in 2015
IEEE 56th Annual Symposium on Foundations of
Computer Science (2015) pp. 792–809.
[22] V. V. Shende, S. S. Bullock, and I. L. Markov,
IEEE Transactions on Computer-Aided Design of
Integrated Circuits and Systems 25, 1000 (2006).
[23] L. K. Grover, Phys. Rev. Lett. 95, 150501 (2005).
[24] A. Ambainis, in 29th International Symposium on
Theoretical Aspects of Computer Science (STACS
2012), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 14, edited by C. Dürr and
T. Wilke (Schloss Dagstuhl–Leibniz-Zentrum fuer
Informatik, Dagstuhl, Germany, 2012) pp. 636–647.
[25] A. Y. Kitaev, arXiv preprint quant-ph/9511026
(1995).
[26] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca,
in Proceedings of the Royal Society of London A:
Mathematical, Physical and Engineering Sciences,
Vol. 454 (The Royal Society, 1998) pp. 339–354.
[27] D. Nagaj, P. Wocjan, and Y. Zhang, Quantum Info.
Comput. 9, 1053 (2009).
[28] G. Wang, arXiv preprint arXiv:1311.1851 (2013).

[29] C. H. Bennett, E. Bernstein, G. Brassard, and
U. Vazirani, SIAM Journal on Computing 26, 1510
(1997).
[30] M. Boyer, G. Brassard, P. Hyer, and A. Tapp,
Fortschritte der Physik 46, 493 (1998).
[31] A. E. Hoerl and R. W. Kennard, Technometrics 12,
55 (1970).
[32] R. Tibshirani, Journal of the Royal Statistical Society. Series B (Methodological) 58, 267 (1996).
[33] G. Sentı́s, J. Calsamiglia, R. Munoz-Tapia, and
E. Bagan, Scientific reports 2, 708 (2012).
[34] H. J. Briegel and G. De las Cuevas, Scientific reports
2, 400 (2012).
[35] G. D. Paparo and M. A. Martin-Delgado, Scientific
Reports 2, 444 (2012).
[36] S. Lloyd, M. Mohseni, and P. Rebentrost, arXiv
preprint arXiv:1307.0411 (2013).
[37] E. Aı̈meur, G. Brassard, and S. Gambs, Machine
Learning 90, 261 (2013).
[38] S. Aaronson, Nature Physics 11, 291 (2015).
[39] J. Adcock, E. Allen, M. Day, S. Frick, J. Hinchliff,
M. Johnson, S. Morley-Short, S. Pallister, A. Price,
and S. Stanisic, arXiv preprint arXiv:1512.02900
(2015).
[40] B. O’Gorman, R. Babbush, A. Perdomo-Ortiz,
A. Aspuru-Guzik, and V. Smelyanskiy, The European Physical Journal Special Topics 224, 163
(2015).
[41] G. D. Paparo, V. Dunjko, A. Makmal, M. A.
Martin-Delgado, and H. J. Briegel, Phys. Rev. X
4, 031002 (2014).
[42] P. Rebentrost, M. Mohseni, and S. Lloyd, Phys.
Rev. Lett. 113, 130503 (2014).
[43] N. Wiebe, C. Granade, C. Ferrie, and D. Cory,
Phys. Rev. A 89, 042314 (2014).
[44] M. Schuld, I. Sinayskiy, and F. Petruccione, Contemporary Physics 56, 172 (2015).
[45] M. Schuld, I. Sinayskiy, and F. Petruccione, Physics
Letters A 379, 660 (2015).
[46] G. Sentı́s, M. Guţă, and G. Adesso, EPJ Quantum
Technology 2, 17 (2015).
[47] Z. Zhao, J. K. Fitzsimons, and J. F. Fitzsimons,
arXiv preprint arXiv:1512.03929 (2015).
[48] S. H. Adachi and M. P. Henderson, arXiv preprint
arXiv:1510.06356 (2015).
[49] N. Wiebe, A. Kapoor, and K. M. Svore, Quantum
Information and Computation 15 (2015).
[50] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, arXiv preprint
arXiv:1611.09347 (2016).
[51] I. Kerenidis and A. Prakash, arXiv preprint
arXiv:1603.08675 (2016).
[52] N. Wiebe, A. Kapoor, and K. M. Svore, Quantum
Information and Computation 16 (2016).
[53] V. Dunjko, J. M. Taylor, and H. J. Briegel, Phys.
Rev. Lett. 117, 130501 (2016).
[54] D. Crawford, A. Levit, N. Ghadermarzy,
J. S. Oberoi, and P. Ronagh, arXiv preprint
arXiv:1612.05695 (2016).

21

[56] M. Benedetti, J. Realpe-Gómez, R. Biswas, and
A. Perdomo-Ortiz, arXiv preprint arXiv:1609.02542
(2016).
[57] M. Benedetti, J. Realpe-Gómez, R. Biswas, and
A. Perdomo-Ortiz, Phys. Rev. A 94, 022308 (2016).

[55] M. H. Amin,
E. Andriyash,
J. Rolfe,
B. Kulchytskyy, and R. Melko, arXiv preprint
arXiv:1601.02036 (2016).

22

