# Quantum Regularized Least Squares

Shantanav Chakraborty  Center for Quantum Science and Technology, IIIT
Hyderabad, Telangana 500032, India  Center for Security, Theory and
Algorithmic Research, IIIT Hyderabad, Telangana 500032, India
<EMAIL_ADDRESS>Aditya Morolia  Center for Quantum Science and Technology,
IIIT Hyderabad, Telangana 500032, India  Center for Computational Natural
Sciences and Bioinformatics, IIIT Hyderabad, Telangana 500032, India
<EMAIL_ADDRESS>Anurudh Peduri  Faculty of Computer
Science, Ruhr University Bochum, 44801 Bochum, Germany  Center for Quantum
Science and Technology, IIIT Hyderabad, Telangana 500032, India  Center for
Security, Theory and Algorithmic Research, IIIT Hyderabad, Telangana 500032,
India<EMAIL_ADDRESS>
###### Abstract

Linear regression is a widely used technique to fit linear models and finds
widespread applications across different areas such as machine learning and
statistics. In most real-world scenarios, however, linear regression problems
are often ill-posed or the underlying model suffers from overfitting, leading
to erroneous or trivial solutions. This is often dealt with by adding extra
constraints, known as regularization. In this paper, we use the frameworks of
block-encoding and quantum singular value transformation (QSVT) to design the
first quantum algorithms for quantum least squares with general
$\ell_{2}$-regularization. These include regularized versions of quantum
ordinary least squares, quantum weighted least squares, and quantum
generalized least squares. Our quantum algorithms substantially improve upon
prior results on quantum ridge regression (polynomial improvement in the
condition number and an exponential improvement in accuracy), which is a
particular case of our result.

To this end, we assume approximate block-encodings of the underlying matrices
as input and use robust QSVT algorithms for various linear algebra operations.
In particular, we develop a variable-time quantum algorithm for matrix
inversion using QSVT, where we use quantum singular value discrimination as a
subroutine instead of gapped phase estimation. This ensures that substantially
fewer ancilla qubits are required for this procedure than prior results. Owing
to the generality of the block-encoding framework, our algorithms are
applicable to a variety of input models and can also be seen as improved and
generalized versions of prior results on standard (non-regularized) quantum
least squares algorithms.

###### Contents

  1. 1 Introduction
    1. 1.1 Linear regression with $\ell_{2}$-regularization
    2. 1.2 Prior work
    3. 1.3 Our contributions
  2. 2 Preliminaries
    1. 2.1 Notation
    2. 2.2 Quantum Input Models
      1. 2.2.1 Quantum Data Structure Input Model
      2. 2.2.2 Sparse Access Input Model
    3. 2.3 Quantum Singular Value Transformation
    4. 2.4 Variable Time Amplitude Amplification
  3. 3 Algorithmic Primitives
    1. 3.1 Arithmetic with Block-Encoded Matrices
    2. 3.2 Robust Quantum Singular Value Discrimination
    3. 3.3 Variable-Time Quantum Algorithm for Matrix Inversion using QSVT
    4. 3.4 Negative Powers of Matrices using QSVT
  4. 4 Quantum Least Squares with General $\ell_{2}$-Regularization
    1. 4.1 Quantum Ordinary Least Squares
    2. 4.2 Quantum Weighted And Generalized Least Squares
      1. 4.2.1 Quantum Weighted Least Squares
      2. 4.2.2 Quantum Generalized Least Squares
  5. 5 Future Directions
  6. A Algorithmic Primitives
    1. A.1 Arithmetic with Block-Encoded Matrices

##  1 Introduction

The problem of fitting a theoretical model to a large set of experimental data
appears across various fields ranging from the natural sciences to machine
learning and statistics [35]. Linear regression is one of the most widely used
procedures for achieving this. By assuming that, for the underlying model,
there exists a linear relationship between a dependent variable and one or
more explanatory variables, linear regression constructs the best linear fit
to the series of data points. Usually, it does so while minimizing the sum of
squared errors - known as the least squares method.

In other words, suppose that we are given $N$ data points
$\\{(a_{i},b_{i})\\}_{i=1}^{N}$ where $\forall
i:a_{i}\in\mathbb{R}^{d},\forall i:b_{i}\in\mathbb{R}$. The assumption is that
each $b_{i}$ is linearly dependent on $a_{i}$ up to some random noise of mean
$0$. Suppose $A$ is the data matrix of dimension $N\times d$, such that its
$i^{\mathrm{th}}$ row is the vector $a_{i}$ and $b\in\mathbb{R}^{N}$ such that
$b=(b_{1},\cdots,b_{N})^{T}$. Then the procedure, known as ordinary least
squares, obtains a vector $x\in\mathbb{R}^{d}$ that minimizes the objective
function $\left|\left|Ax-b\right|\right|^{2}_{2}$. This problem has a closed-
form solution given by $x=(A^{T}A)^{-1}A^{T}b=A^{+}b$, where $A^{+}$ denotes
the Moore-Penrose inverse of the matrix $A$. Thus computationally, finding the
best fit by linear regression reduces to finding the pseudoinverse of a matrix
that represents the data, a task that is expensive for classical machines for
large data sets.

In practice, however, least squares regression runs into problems such as
overfitting. For instance, the solution might fit most data points, even those
corresponding to random noise. Furthermore, the linear regression problem may
also be ill-posed, for instance, when the number of variables exceeds the
number of data points rendering it impossible to fit the data. These issues
come up frequently with linear regression models and result in erroneous or
trivial solutions. Furthermore, another frequent occurrence is that the data
matrix $A$ has linearly dependent columns. In this scenario, the matrix
$A^{T}A$ is not full rank and therefore is not invertible.

Regularization is a widely used technique to remedy these problems, not just
for linear regression but for inverse problems, in general [11]. In the
context of linear regression, broadly, this involves adding a penalty term to
the objective function, which constrains the solution of the regression
problem. For instance, in the case of $\ell_{2}$-regularization, the objective
is to obtain $x$ that minimizes

$\norm{Ax-b}^{2}_{2}+\lambda\norm{Lx}^{2}_{2}$ (1)

where $L$ is an appropriately chosen penalty matrix (or regularization matrix)
of dimension $N\times d$ and $\lambda>0$ is the regularization parameter, an
appropriately chosen constant. This regularization technique is known as
general $\mathit{\ell_{2}}$-regularization or Tikhonov regularization in the
literature [17, 18, 4, 12, 43]. It is a generalization of ridge regression
which corresponds to the case when $L$ is the identity matrix [20, 33, 42].
The closed-form solution of the general $\ell_{2}$-regularized ordinary least
squares problem is given by

$x=\left(A^{T}A+\lambda L^{T}L\right)^{-1}A^{T}b.$ (2)

A straightforward observation is that even when $A^{T}A$ is singular, a
judicious choice of the penalty matrix $L$ can ensure that the effective
condition number (ratio of the maximum and the minimum non-zero singular
values) of the overall matrix is finite and $A^{T}A+\lambda L^{T}L$ is
invertible.

In this paper, we develop quantum algorithms for linear regression with
general $\ell_{2}$-regularization. If the optimal solution is
$x=(x_{1},\cdots,x_{d})^{T}$, then our quantum algorithm outputs a quantum
state that is $\delta$-close to $\ket{x}=\sum_{j=1}^{d}x_{j}\ket{j}/\norm{x}$,
assuming access to the matrices $A,L$, and the quantum state $\ket{b}$ via
general quantum input models.

In several practical scenarios, depending on the underlying theoretical model,
generalizations of the ordinary least squares (OLS) technique are more useful
to fit the data. For instance, certain samples may be of more importance (and
therefore have more weight) than the others, in which case weighted least
squares (WLS) is preferred. Generalized least squares (GLS) is used when the
underlying samples obtained are correlated. These techniques also suffer from
the issues commonplace with OLS, warranting the need for regularization [43].
Consequently, we also design algorithms for quantum WLS with general
$\ell_{2}$-regularization and quantum GLS with general
$\ell_{2}$-regularization.

Organization of the paper: In the remainder of Section 1, we formally describe
$\ell_{2}$-regularized versions of OLS, WLS, and GLS (Section 1.1), discuss
prior and related work (Section 1.2), and outline our contributions and
results (Section 1.3). In Section 2, we briefly outline the framework of
block-encoding and quantum input models that are particular instances of it
(Section 2.2). We also briefly introduce quantum singular value transformation
(QSVT) (Section 2.3) and variable time amplitude amplification (VTAA) (Section
2.4). Following this, in Section 3, we develop several algorithmic primitives
involving arithmetic of block-encodings (Section 3.1), quantum singular value
discrimination (Section 3.2) and quantum linear algebra using QSVT (Section
3.3). These are the technical building blocks for designing our quantum
regularized regression algorithms. Using these algorithmic primitives, we
design quantum algorithms for the quantum least squares with
$\ell_{2}$-regularization in Section 4. Finally, we conclude by discussing
some possible future research directions in Section 5.

###  1.1 Linear regression with $\ell_{2}$-regularization

Suppose we are given data points $\\{(a_{i},b_{i})\\}_{i=1}^{N}$, where
$\forall i:a_{i}\in\mathbb{R}^{d},\forall i:b_{i}\in\mathbb{R}$ such that
$(a_{i},b_{i})\sim_{i.i.d}\mathcal{D}$, i.e. they are sampled i.i.d. from some
unknown distribution $\mathcal{D}$, assumed to be linear. We want to find a
vector $x\in\mathbb{R}^{d}$ such that the inner product $x^{T}a_{j}$ is a good
predictor for the target $b_{j}$ for some unknown $a_{j}$. This can be done by
minimizing the total squared loss over the given data points,

$\mathcal{L}_{O}:=\sum_{j}(x^{T}a_{j}-b_{j})^{2},$ (3)

leading to the ordinary least squares (OLS) optimization problem. The task
then is to find $x\in\mathbb{R}^{d}$ that minimizes $\norm{Ax-b}_{2}^{2}$,
where $A$ is the $N\times d$ data matrix such that the $i^{th}$ row of $A$ is
$a_{i}$, and the $i^{th}$ element of the vector $b$ is $b_{i}$. Assuming that
$A^{T}A$ is non-singular, the optimal $x$ satisfies

$x=(A^{T}A)^{-1}A^{T}b=A^{+}b,$ (4)

which corresponds to solving a linear system of equations.

Suppose that out of the samples present in the data, we have higher confidence
in some of them than others. In such a scenario, the $i^{\mathrm{th}}$
observation can be assigned a weight $w_{i}\in\mathbb{R}$. This leads to a
generalization of the OLS problem to weighted least squares (WLS). In order to
obtain the best linear fit, the task is now to minimize the weighted version
of the loss

$\mathcal{L}_{W}:=\sum_{j}w_{j}(x^{T}a_{j}-b_{j})^{2}.$ (5)

As before, assuming $A^{T}WA$ is non-singular, the above loss function has the
following closed-form solution:

$x=(A^{T}WA)^{-1}A^{T}Wb,$ (6)

where $W$ is a diagonal matrix with $w_{i}$ being the $i^{\mathrm{th}}$
diagonal element.

There can arise scenarios where there exists some correlation between any two
samples. For generalized least squares (GLS), the presumed correlations
between pairs of samples are given in a symmetric, non-singular covariance
matrix $\Omega$. This objective is to find the vector $x$ that minimizes

$\mathcal{L}_{\Omega}:=\sum_{i,j}(\Omega^{-1})_{i,j}(x^{T}a_{i}-b_{i})(x^{T}a_{j}-b_{j}).$
(7)

Similarly, the closed-form solution for GLS is given by

$x=(A^{T}\Omega^{-1}A)^{-1}A^{T}\Omega^{-1}b.$ (8)

As mentioned previously, in several practical scenarios, the linear regression
problem may be ill-posed or suffer from overfitting. Furthermore, the data may
be such that some of the columns of the matrix $A$ are linearly dependent.
This shrinks the rank of $A$, and consequently of the matrix $A^{T}A$,
rendering it singular and, therefore non-invertible. Recall that the closed-
form solution of OLS exists only if $A^{T}A$ is non-singular, which is no
longer the case. Such scenarios arise even for WLS and GLS problems [43].

In such cases, one resorts to regularization to deal with them. Let
$\mathcal{L}$ be the loss function to be minimized for the underlying least
squares problem (such as OLS, WLS, or GLS). Then general
$\ell_{2}$-regularization (Tikhonov regularization) involves an additional
penalty term so that the objective now is to find the vector
$x\in\mathbb{R}^{d}$ that minimizes

$\mathcal{L}+\lambda\norm{Lx}^{2}_{2}.$ (9)

Here $\lambda$, known as the regularization parameter, is a positive constant
that controls the size of the vector $x$, while $L$ is known as the penalty
matrix (or regularization matrix) that defines a (semi)norm on the solution
through which the size is measured. The solution to the Tikhonov
regularization problem also has a closed-form solution. For example, in the
OLS problem, when $\mathcal{L}=\mathcal{L}_{O}$, we have that

$x=(A^{T}A+\lambda L^{T}L)^{-1}A^{T}b.$ (10)

It is worth noting that when $L=I$, the $\ell_{2}$-regularized OLS problem is
known as ridge regression. For the unregularized OLS problem, the singular
values of $A$, $\sigma_{j}$ are mapped to $1/\sigma_{j}$. The penalty term due
to $\ell_{2}$-regularization, results in a shrinkage of the singular values.
This implies that even in the scenario where $A$ has linearly dependent
columns (some $\sigma_{j}=0$) and $(A^{T}A)^{-1}$ does not exist, the inverse
$(A^{T}A+\lambda L^{T}L)^{-1}$ is well defined for $\lambda>0$ and any
positive-definite $L$. Throughout this article, we refer to such an $L$ (which
is positive definite) as a good regularizer. The penalty matrix $L$ allows for
penalizing each regression parameter differently and leads to joint shrinkage
among the elements of $x$. It also determines the rate and direction of
shrinkage. In the special case of ridge regression, as $L=I$, the penalty
shrinks each element of $x$ equally along the unit vectors $e_{j}$. Also note
that by definition, $I$ is a good regularizer.

Closed-form expressions can also be obtained for the WLS and the GLS problem
($\mathcal{L}=\mathcal{L}_{W},\mathcal{L}_{\Omega}$ respectively), and finding
the optimal solution $x$ reduces to solving a linear system. The quantum
version of these algorithms output a quantum state that is $\epsilon$-close
$\ket{x}=\sum_{j}x_{j}\ket{j}/\norm{x}$.

Throughout this work, while designing our quantum algorithms, we shall assume
access (via a block-encoding) to the matrices $A$, $W$, $\Omega$, and $L$ and
knowledge of the parameter $\lambda$. Classically, however, the regularization
matrix $L$ and the optimal parameter $\lambda$ are obtained via several
heuristic techniques [18, 12, 43].

###  1.2 Prior work

Quantum algorithms for (unregularized) linear regression was first developed
by Wiebe et al. [45], wherein the authors made use of the HHL algorithm for
solving a linear system of equations [19]. Their algorithm assumes query
access to a sparse matrix $A$ (sparse-access-model) and to a procedure to
prepare $\ket{b}=\sum_{i}b_{i}\ket{i}$. They first prepare a quantum state
proportional to $A^{T}\ket{b}$, and then use the HHL algorithm to apply the
operator $(A^{T}A)^{-1}$ to it. Overall the algorithm runs in a time scaling
as $\kappa^{6}_{A}$ (the condition number of $A$) and inverse polynomial in
the accuracy $\delta$. Subsequent results have considered the problem of
obtaining classical outputs for linear regression. For instance, in Ref. [44],
$A^{+}$ is directly applied to the quantum state $\ket{b}$, followed by
amplitude estimation to obtain the entries of $x$. On the other hand, Ref.
[39] used the techniques of quantum principal component analysis in [27] to
predict a new data point for the regression problem. These algorithms also
work in the sparse access model and run in a time that scales as
$\mathrm{poly}\left(\kappa,1/\delta\right)$. Kerenidis and Prakash [23]
provided a quantum algorithm for the WLS problem wherein they used a classical
data structure to store the entries of $A$ and $W$. Furthermore, they assumed
QRAM access to this data structure [36, 22] that would allow the preparation
of quantum states proportional to the entries of $A$ and $W$ efficiently. They
showed that in this input model (quantum data structure model), an iterative
quantum linear systems algorithm can prepare $\ket{x}$ in time
$\widetilde{O}(\mu\kappa^{3}/\delta)$, where $\kappa$ is the condition number
of the matrix $A^{T}\sqrt{W}$ while $\mu=\norm{\sqrt{W}A}_{F}$. Chakraborty et
al. [7] applied the framework of block-encoding along with (controlled)
Hamiltonian simulation of Low and Chuang [26] to design improved quantum
algorithms for solving linear systems. Quantum algorithms developed in the
block-encoding framework are applicable to a wide variety of input models,
including the sparse access model and the quantum data structure model of
[23]. They applied their quantum linear systems solver to develop quantum
algorithms for quantum weighted least squares and generalized least squares.
Their quantum algorithm for WLS has a complexity that is in
$\widetilde{O}\left(\alpha\kappa\mathrm{polylog}(Nd/\delta)\right)$, where
$\alpha=s$, the sparsity of the matrix $A^{T}\sqrt{W}$ in the sparse access
model while $\alpha=\norm{\sqrt{W}A}_{F}$, for the quantum data structure
input model. For GLS, their quantum algorithm outputs $\ket{x}$ in cost
$\widetilde{O}\left(\kappa_{A}\kappa_{\Omega}(\alpha_{A}+\alpha_{\Omega}\kappa_{\Omega})\mathrm{polylog}(1/\delta)\right)$,
where $\kappa_{A}$ and $\kappa_{\Omega}$ are the condition numbers of $A$ and
$\Omega$ respectively while $\alpha_{A}$ and $\alpha_{\Omega}$ are parameters
that depend on how the matrices $A$ and $\Omega$ are accessed in the
underlying input model.

While quantum linear regression algorithms have been designed and subsequently
improved over the years, quantum algorithms for regularized least squares have
not been developed extensively. Yu et al. [46] developed a quantum algorithm
for ridge regression in the sparse access model using the LMR scheme [27] for
Hamiltonian simulation and quantum phase estimation, which they then used to
determine the optimal value of the parameter $\lambda$. Their algorithm to
output $\ket{x}$ has a cubic dependence on both $\kappa$ and $1/\delta$. They
use this as a subroutine to determine a good value of $\lambda$. A few other
works [40, 10] have considered the quantum ridge regression problem in the
sparse access model, all of which can be implemented with
$\mathrm{poly}(\kappa,1/\delta)$ cost.

Recently, Chen and de Wolf designed quantum algorithms for lasso
($\ell_{1}$-regularization) and ridge regressions from the perspective of
empirical loss minimization [5]. For both lasso and ridge, their quantum
algorithms output a classical vector $\widetilde{x}$ whose loss (mean squared
error) is $\delta$-close to the minimum achievable loss. In this context, they
prove a quantum lower bound of $\Omega(d/\delta)$ for ridge regression which
indicates that in their setting, the dependence on $d$ cannot be improved on a
quantum computer (the classical lower bound is also linear in $d$ and there
exists a matching upper bound). Note that $\widetilde{x}$ is not necessarily
close to the optimal solution $x$ of the corresponding least squares problem,
even though their respective loss values are. Moreover, their result (of
outputting a classical vector $\widetilde{x}$) is incomparable to our
objective of obtaining a quantum state encoding the optimal solution to the
regularized regression problem.

Finally, Gilyén et al. obtained a “dequantized” classical algorithm for ridge
regression assuming norm squared access to input data similar to the quantum
data structure input model [15]. Furthermore, similar to the quantum setting
where the output is the quantum state $\ket{x}=\sum_{j}x_{j}\ket{j}/\norm{x}$
instead of $x$ itself, their algorithm obtains samples from the distribution
$x_{j}^{2}/\norm{x}^{2}$. For the regularization parameter
$\lambda=O\left(\norm{A}\norm{A}_{F}\right)$, the running time of their
algorithm is in $\widetilde{O}\left(\kappa^{12}r_{A}^{3}/\delta^{4}\right)$,
where $r_{A}$ is the rank of $A$. Their result (and several prior results)
does not have a polynomial dependence on the dimension of $A$ and therefore
rules out the possibility of generic exponential quantum speedup (except in
$\delta$) in the quantum data structure input model.

###  1.3 Our contributions

In this work, we design the first quantum algorithms for OLS, WLS, and GLS
with general $\ell_{2}$-regularization. We use the Quantum Singular Value
Transformation (QSVT) framework introduced by Gilyén et al [14]. We assume
that the relevant matrices are provided as input in the block-encoding model,
in which access to an input matrix $A$ is given by a unitary $U_{A}$ whose
top-left block is (close to) $A/\alpha$. The parameter $\alpha$ takes specific
values depending on the underlying input model. QSVT then allows us to
implement nearly arbitrary polynomial transformations to a block of a unitary
matrix using a series of parameterized, projector-controlled rotations
(quantum signal processing [25]).

More precisely, given approximate block-encodings of the data matrix $A$ and
the regularizing matrix $L$, and a unitary procedure to prepare the state
$\ket{b}$, our quantum algorithms output a quantum state that is
$\delta$-close to $\ket{x}$, the quantum state proportional to the
$\ell_{2}$-regularized ordinary least squares (or weighted least squares or
generalized least squares problem). We briefly summarize the query
complexities of our results in Table 1.

For the OLS problem with general $\ell_{2}$-regularization (Section 4.2,
Theorem 32), we design a quantum algorithm which given an
$(\alpha_{A},a_{A},\varepsilon_{A})$-block-encoding of $A$ (implemented in
cost $T_{A}$), an $(\alpha_{L},a_{L},\varepsilon_{L})$-block-encoding of $L$
(implemented in cost $T_{L}$), a parameter $\lambda>0$, and a procedure to
prepare $\ket{b}$ (in cost $T_{b}$), outputs a quantum state which is
$\delta$-close to $\ket{x}$. The algorithm has a cost

$\order{\kappa\log\kappa\left(\left(\frac{\alpha_{A}+\sqrt{\lambda}\alpha_{L}}{\norm{A}+\sqrt{\lambda}\norm{L}}\right)\log\left(\frac{\kappa}{\delta}\right)\left(T_{A}+T_{L}\right)+T_{b}\right)}$

where $\kappa$ can be thought of as a modified condition number, related to
the effective condition numbers of $A$ and $L$. When $L$ is a good
regularizer, this is given by the expression

$\kappa=\kappa_{L}\left(1+\frac{\norm{A}}{\sqrt{\lambda}\norm{L}}\right),$

Notice that $\kappa$ is independent of $\kappa_{A}$, the condition number of
the data matrix $A$, which underscores the advantage of regularization. The
parameters $\alpha_{A}$ and $\alpha_{L}$ take specific values depending on the
underlying input model. For the sparse access input model, $\alpha_{A}=s_{A}$
and $\alpha_{L}=s_{L}$, the respective sparsities of the matrices $A$ and $L$.
On the other hand for the quantum data structure input model,
$\alpha_{A}=\norm{A}_{F}$ and $\alpha_{L}=\norm{L}_{F}$. Consequently, the
complexity of Quantum Ridge Regression can be obtained by substituting $L=I$
in the above complexity as

$\order{\log\kappa\left(\frac{\alpha_{A}}{\sqrt{\lambda}}\log\left(\frac{\kappa}{\delta}\right)T_{A}+\kappa
T_{b}\right)}$

where $\kappa=1+\norm{A}/\sqrt{\lambda}$, by noting that the block-encoding of
$I$ is trivial while the norm and condition number of the identity matrix is
one. For this problem of quantum ridge regression, our quantum algorithms are
substantially better than prior results [40, 46, 10], exhibiting a polynomial
improvement in $\kappa$ and an exponential improvement in $1/\delta$.

For the $\ell_{2}$-regularized GLS problem (Section 4, Theorem 42), we design
a quantum algorithm that along with approximate block-encodings of $A$ and
$L$, takes as input an
$(\alpha_{\Omega},a_{\Omega},\varepsilon_{\Omega})$-lock-encoding of the
matrix $\Omega$ (implementable at a cost of $T_{\Omega}$) to output a state
$\delta$-close to $\ket{x}$ at a cost of

$\order{\kappa\sqrt{{\kappa_{\Omega}}}\log\kappa\left(\left(\frac{\alpha_{A}}{\norm{A}}T_{A}+\frac{\alpha_{L}}{\norm{L}}T_{L}+\frac{\alpha_{\Omega}{\kappa_{\Omega}}}{{\norm{\Omega}}}T_{\Omega}\right)\log^{3}\left(\frac{\kappa{\kappa_{\Omega}}\norm{A}\norm{L}}{\delta{\norm{\Omega}}}\right)+T_{b}\right)}$

In the above complexity, when $L$ is a good regularizer, the modified
condition number $\kappa$ is defined as

$\kappa=\kappa_{L}\left(1+\frac{\sqrt{{\kappa_{\Omega}}}\norm{A}}{\sqrt{\lambda\norm{\Omega}}\norm{L}}\right)$

The WLS problem is a particular case of GLS, wherein the matrix $\Omega$ is
diagonal. However, we show that better complexities for the
$\ell_{2}$-regularized WLS problem can be obtained if we assume QRAM access to
the diagonal entries of $W$ (Section 4, Theorem 39 and Theorem 40).

Table 1 summarizes the complexities of our algorithms for quantum linear
regression with general $\ell_{2}$-regularization. For better exposition, here
we assume that $\norm{A},\norm{L},\norm{\Omega}$ and
$\lambda=\Theta\left(1\right)$. For the general expression of the
complexities, we refer the readers to Section 4.

Problem | Unregularized |  $\mathbf{\ell_{2}}$-Regularized
---|---|---
Quantum OLS | $\widetilde{\mathcal{O}}\left(\alpha_{A}\kappa_{A}\log\left(1/\delta\right)\right)$ | $\widetilde{\mathcal{O}}\left((\alpha_{A}+\alpha_{L})\kappa_{L}\log\left(1/\delta\right)\right)$
Quantum GLS | $\widetilde{\mathcal{O}}\left(\left(\alpha_{A}+\alpha_{\Omega}\kappa_{\Omega}\right)\kappa_{A}\sqrt{\kappa_{\Omega}}\log^{3}\left(1/\delta\right)\right)$ | $\widetilde{\mathcal{O}}\left(\left(\alpha_{A}+\alpha_{L}+\alpha_{\Omega}\kappa_{\Omega}\right)\kappa_{L}\sqrt{\kappa_{\Omega}}\log^{3}\left(1/\delta\right)\right)$
Table 1: Complexity of quantum linear regression algorithms with and without
general $\ell_{2}$-regularization. All of these algorithms require only
$\Theta({\log\kappa})$ additional qubits.

In order to derive our results, we take advantage of the ability to
efficiently perform arithmetic operations on block-encoded matrices, as
outlined in Section 3. Along with this, we use QSVT to perform linear
algebraic operations on block-encoded matrices. To this end, adapt the results
in Refs. [14, 34] to our setting. One of our contributions is that we work
with robust versions of many of these algorithms. In prior works, QSVT is
often applied to block-encoded matrices, assuming perfect block-encoding. For
the quantum algorithms in this paper, we rigorously obtain the precision
$\varepsilon$ required to obtain a $\delta$-approximation of the desired
output state.

For instance, a key ingredient of our algorithm for regularized least squares
is to make use of QSVT to obtain $A^{+}$, given an $\varepsilon$-approximate
block-encoding of $A$. In order to obtain a (near) optimal dependence on the
condition number of $A$ by applying variable-time amplitude amplification
(VTAA) [3], we recast the standard QSVT algorithm as a variable stopping-time
quantum algorithm. Using QSVT instead of controlled Hamiltonian simulation
ensures that the variable-time quantum procedure to prepare $A^{+}\ket{x}$ has
a slightly better running time (by a log factor) and considerably fewer
additional qubits than Refs. [9, 7].

Furthermore, for the variable time matrix inversion algorithm, a crucial
requirement is the application of the inversion procedure to the portion of
the input state that is spanned by singular values larger than a certain
threshold. In order to achieve this, prior results have made use of Gapped
Phase Estimation (GPE), which is a simple variant of the standard phase
estimation procedure that decides whether the eigenvalue of a Hermitian
operator is above or below a certain threshold [3, 9, 7]. However, GPE can
only be applied to a Hermitian matrix and requires additional registers that
store the estimates of the phases, which are never used for variable-time
amplitude amplification. In this work, instead of GPE, we develop a robust
version of quantum singular value discrimination (QSVD) using QSVT, which can
be directly applied to non-Hermitian matrices. This algorithm decides whether
some singular value of a matrix is above or below a certain threshold without
storing estimates of the singular values. This leads to a space-efficient
variable time quantum algorithm for matrix inversion by further reducing the
number of additional qubits required by a factor of
$O(\log^{2}(\kappa/\delta))$ as compared to prior results [9, 7].
Consequently, this also implies that in our framework, quantum algorithms for
(unregularized) least squares (which are special cases of our result) have
better complexities than those of Ref. [7].

##  2 Preliminaries

This section lays down the notation, and introduces the quantum singular value
transformation (QSVT) and block-encoding frameworks, which are used to design
the algorithm for quantum regression.

###  2.1 Notation

For a matrix $A\in\mathbb{R}^{N\times d}$, $A_{i,.}$ denotes the $i^{th}$ row
of $A$, and $\norm{A_{i,\cdot}}$ denotes the vector norm of $A_{i,.}^{T}$.
$s^{A}_{r}$ and $s^{A}_{c}$ denote the row and column sparsity of the matrix,
which is the maximum number of non-zero entries in any row and any column,
respectively.

Singular Value Decomposition. The decomposition $A=W\Sigma V^{\dagger}$, where
$W$ and $V$ are unitary and $\Sigma$ is a diagonal matrix, represents the
singular value decomposition (SVD) of $A$. All matrices can be decomposed in
this form. The diagonal entries of $\Sigma$, usually denoted by
$\sigma(A)=\\{\sigma_{j}\\}$, is the multiset of all singular values of $A$,
which are real and non-negative. $\sigma_{\max}$ and $\sigma_{\min}$ denote
the maximum and minimum singular values of $A$. $r(A)=\mathrm{rank}(A)$ is the
number of non-zero singular values of $A$. The columns of
$W,\leavevmode\nobreak\ V$ (denoted by $\\{\ket{w_{j}}\\}$ and
$\\{\ket{v_{j}}\\}$) are the left and right singular vectors of $A$. Thus
$A=\sum_{j}^{r}\sigma_{j}\ket{w_{j}}\bra{v_{j}}$. The singular vectors of $A$
can be computed as the positive square roots of the eigenvalues of
$A^{\dagger}A$ (which is positive semi-definite and therefore has non-negative
real eigenvalues.)

Effective Condition Number. $\kappa_{A}$ denotes (an upper bound on) the
effective condition number of $A$, defined as the ratio of the maximum and
minimum non-zero singular values of $A$. Let $\sigma_{\max}\left(A\right)$ be
the largest singular value of $A$, and $\sigma_{\min}\left(A\right)$ be the
smallest singular value of $A$. Additionally, let
$\widetilde{\sigma}_{\min}\left(A\right)$ be the smallest non-zero singular
value of $A$. Then

$\kappa_{A}\geq\frac{\sigma_{\max}\left(A\right)}{\widetilde{\sigma}_{\min}\left(A\right)}=\sqrt{\frac{\lambda_{\max}(A^{\dagger}A)}{\widetilde{\lambda}_{\min}(A^{\dagger}A)}}$

If $A$ is full-rank, then
$\widetilde{\sigma}_{\min}\left(A\right)=\sigma_{\min}\left(A\right)$, and
$\kappa_{A}$ becomes the condition number of the matrix. In this text, unless
stated otherwise, we always refer to $\kappa_{A}$ as (an upper bound on)
effective condition number of a matrix, and not the true condition number.

Norm. Unless otherwise specified, $\norm{A}$ denotes the spectral norm of $A$,
while $\norm{A}_{F}$ denotes the Frobenius norm of $A$, defined as

$\displaystyle\norm{A}:=\max_{x\neq
0}\frac{\norm{Ax}}{\norm{x}}=\sigma_{\max}(A)$
$\displaystyle\norm{A}_{F}:=\sqrt{\sum_{j=1}^{r}\sigma_{j}^{2}}$

Unless otherwise specified, when $A$ is assumed to be normalized, it is with
respect to the spectral norm.

Soft-O Complexity. Finally, we use $f=\widetilde{\mathcal{O}}\left(g\right)$
to denote $f=\order{g\cdot\mathrm{polylog}(g)}$.

Controlled Unitaries. If $U$ is a $s$-qubit unitary, then $C\text{-}U$ is a
$(s+1)$-qubit unitary defined by

$C\text{-}U=\outerproduct{0}{0}\otimes I_{s}+\outerproduct{1}{1}\otimes U$

Throughout this text whenever we state that the time taken to implement a
unitary $U_{A}$ is $T_{A}$ and the cost of an algorithm is $\order{nT_{A}}$,
we imply that the algorithm makes $n$ uses of the unitary $U_{A}$. Thus, if
the circuit depth of $U_{A}$ is $T_{A}$, the circuit depth of our algorithm is
$\order{nT_{A}}$.

###  2.2 Quantum Input Models

The complexities of quantum algorithms often depend on how the input data is
accessed. For instance, in quantum algorithms for linear algebra (involving
matrix operations), it is often assumed that there exists a black-box that
returns the positions of the non-zero entries of the underlying matrix when
queried. The algorithmic running time is expressed in terms of the number of
queries made to this black-box. Such an input model, known as the Sparse
Access Model, helps design efficient quantum algorithms whenever the
underlying matrices are sparse. Various other input models exist, and quantum
algorithms are typically designed and optimized for specific input models.

Kerenidis and Prakash [22, Section 5.1] introduced a different input model,
known as the quantum data structure model, which is more conducive for
designing quantum machine learning algorithms. In this model, the input data
(e.g: entries of matrices) arrive online and are stored in a classical data
structure (often referred to as the KP-tree in the literature), which can be
queried in superposition by using a QRAM. This facilitates efficiently
preparing quantum states corresponding to the rows of the underlying matrix,
that can then be used for performing several matrix operations. Subsequently,
several quantum-inspired classical algorithms have also been developed
following the breakthrough result of Tang [41]. Such classical algorithms have
the same underlying assumptions as the quantum algorithms designed in the data
structure input model and are only polynomially slower provided the underlying
matrix is low rank.

In this work, we will consider the framework of block-encoding, wherein it is
assumed that the input matrix $A$ (up to some sub-normalization) is stored in
the left block of some unitary. The advantage of the block-encoding framework,
which was introduced in a series of works [26, Definition 1], [7, Section 1],
[14, Section 1.3], is that it can be applied to a wide variety of input
models. For instance, it can be shown that both the sparse access input model
as well as the quantum data structure input model are specific instances of
block-encoded matrices [7, Sections 2.2 and 2.4], [14, Section 5.2]. Here we
formally define the framework of block-encoding and also express the sparse
access model as well as the quantum data structure model as block-encodings.
We refer the reader to [7, 14] for proofs.

######  Definition 1 (Block Encoding, restated from [14], Definition 24).

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

$\norm{A-\alpha(\bra{0}^{\otimes a}\otimes I)U_{A}(\ket{0}^{\otimes a}\otimes
I)}\leq\varepsilon.$ (11)

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

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

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

In the subsequent sections, we provide an outline of the quantum data
structure model and the sparse access model which are particular instances of
the block encoding framework.

####  2.2.1 Quantum Data Structure Input Model

Kerenidis and Prakash introduced a quantum accessible classical data structure
which has proven to be quite useful for designing several quantum algorithms
for linear algebra [22]. The classical data structure stores entries of
matrices or vectors and can be queried in superposition using a QRAM (quantum
random access memory). We directly state the following theorem from therein.

######  Theorem 2 (Implementing quantum operators using an efficient data
structure, [36, 22]).

Let $A\in\mathbb{R}^{N\times d}$, and $w$ be the number of non-zero entries of
$A$. Then there exists a data structure of size
$\order{w\log^{2}\left(dN\right)}$ that given the matrix elements
$(i,j,a_{ij})$, stores them at a cost of $\order{\log\left(dN\right)}$
operations per element. Once all the non-zero entries of $A$ have been stored
in the data structure, there exist quantum algorithms that are
$\varepsilon$-approximations to the following maps:

$U:\ket{i}\ket{0}\mapsto\frac{1}{\norm{A_{i,\cdot}}}\sum_{j=1}^{d}a_{i,j}\ket{i,j}=\ket{\psi_{i}},$
$V:\ket{0}\ket{j}\mapsto\frac{1}{\norm{A}_{F}}\sum_{i=1}^{N}\norm{A_{i},.}\ket{i,j}=\ket{\phi_{j}}$

where $\norm{A_{i,\cdot}}$ is the norm of the $i^{\mathrm{th}}$ row of $A$ and
the second register of $\ket{\psi_{i}}$ is the quantum state corresponding to
the $i^{\mathrm{th}}$ row of $A$. These operations can be applied at a cost of
$\order{\mathrm{polylog}(Nd/\varepsilon)}$.

It was identified in Ref. [7] that if a matrix $A$ is stored in this quantum
accessible data structure, there exists an efficiently implementable block-
encoding of $A$. We restate their result here.

######  Lemma 3 (Implementing block encodings from quantum data structures,
[7], Theorem 4).

Let the entries of the matrix $A\in\mathbb{R}^{N\times d}$ be stored in a
quantum accessible data structure, then there exist unitaries $U_{R},U_{L}$
that can be implemented at a cost of
$\order{\mathrm{polylog}(dN/\varepsilon)}$ such that $U_{R}^{\dagger}U_{L}$ is
a $(\norm{A}_{F},\lceil\log\left(d+N\right)\rceil,\varepsilon)$-block-encoding
of $A$.

###### Proof.

The unitaries $U_{R}$ and $U_{L}$ can be implemented via $U$ and $V$ in the
previous lemma. Let $U_{R}=U$ and $U_{L}=V.\texttt{SWAP}$. Then for
$s=\lceil\log(d+N)\rceil$ we have

$U_{R}:\leavevmode\nobreak\ \ket{i}\ket{0^{s}}\rightarrow\ket{\psi_{i}},$

and

$U_{L}:\leavevmode\nobreak\ \ket{j}\ket{0^{s}}\rightarrow\ket{\phi_{j}},$

So we have that the top left block of $U^{{\dagger}}_{R}U_{L}$ is

$\sum_{i=1}^{N}\sum_{j=1}^{d}\braket{\psi_{i}}{\phi_{j}}\ket{i,0}\bra{j,0}$

Now

$\displaystyle\braket{\psi_{i}}{\phi_{j}}$
$\displaystyle=\sum_{k=1}^{d}\sum_{\ell=1}^{N}\dfrac{a_{ik}}{\norm{A_{i,\cdot}}}\cdot\dfrac{\norm{A_{\ell}}}{\norm{A}_{F}}\underbrace{\braket{i,k}{l,j}}_{:=\delta_{i,l}.\delta_{k,j}}$
$\displaystyle=\dfrac{a_{ij}}{\norm{A}_{F}}.$

Moreover since only $\varepsilon$-approximations of $U$ and $V$ can be
implemented we have that $U^{{\dagger}}_{R}U_{L}$ is a
$(\norm{A}_{F},\lceil\log(n+d)\rceil,\varepsilon)$ block encoding of $A$
implementable with the same cost as $U$ and $V$. ∎

In Ref. [23] argued that in certain scenarios, storing the entries of
$A^{(p)},(A^{1-p})^{\dagger}$ might be useful as compared to storing $A$, for
some $p\in[0,1]$. In such cases, the quantum data structure is a
$(\mu_{p},\lceil\log(N+d)\rceil,\varepsilon)$ block encoding of $A$, where
$\mu_{p}(A)=\sqrt{s_{2p}(A).s_{2(1-p)}(A^{T})}$ such that
$s_{p}(A):=\max_{j}\norm{A_{j,\cdot}}_{q}^{q}$. Throughout the work, whenever
our results are expressed in the quantum data structure input model, we shall
state our complexity in terms of $\mu_{A}$. When the entries of $A$ are
directly stored in the data structure, $\mu_{A}=\norm{A}_{F}$. Although, we
will not state it explicitly each time, our results also hold when fractional
powers of $A$ are stored in the database and simply substituting
$\mu_{A}=\mu_{p}(A)$, yields the required complexity.

####  2.2.2 Sparse Access Input Model

The sparse access input model considers that the input matrix
$A\in\mathbb{R}^{N\times d}$ has row sparsity $s_{r}$ and column sparsity
$s_{c}$. Furthermore, it assumes that the entries of $A$ can be queried via an
oracle as

$O_{A}:\ket{i}\ket{j}\ket{0}^{\otimes
b}\mapsto\ket{i}\ket{j}\ket{a_{ij}}\quad\forall i\in[N],j\in[d],$

and the indices of the non-zero elements of each row and column can be queried
via the following oracles:

$\displaystyle O_{r}:\ket{i}\ket{j}\mapsto\ket{i}\ket{r_{ij}}\quad\forall
i\in[N],k\in[s_{r}],$ $\displaystyle
O_{c}:\ket{i}\ket{j}\mapsto\ket{c_{ij}}\ket{j}\quad\forall
i\in[d],k\in[s_{c}]$

where $r_{ij}$ is the $j^{\mathrm{th}}$ non-zero entry of the
$i^{\mathrm{th}}$ row of $A$ and $c_{ij}$ is the $i^{\mathrm{th}}$ non-zero
entry of the $j^{\mathrm{th}}$ column of $A$. Gilyén et al. [14] showed that a
block encoding of a sparse $A$ can be efficiently prepared by using these
three oracles. We restate their lemma below.

######  Lemma 4 (Constructing a block-encoding from sparse-access to matrices,
[13], Lemma 48).

Let $A\in\mathbb{R}^{N\times d}$ be an $s_{r},s_{c}$ row, column sparse matrix
given as a sparse access input. Then for all $\varepsilon\in(0,1)$, we can
implement a
$(\sqrt{s_{c}s_{r}},\mathrm{polylog}(Nd/\varepsilon),\varepsilon)$-block-
encoding of $A$ with $\order{1}$ queries to $O_{r},O_{c},O_{A}$ and
$\mathrm{polylog}(Nd/\varepsilon)$ elementary quantum gates.

Throughout the paper, we shall assume input matrices are accessible via
approximate block-encodings. This also allows us to write down the
complexities of our quantum algorithms in this general framework.
Additionally, we state the complexities in both the sparse access input model
as well as the quantum accessible data structure input model as particular
cases.

###  2.3 Quantum Singular Value Transformation

In a seminal work, Gilyén et al. presented a framework to apply an arbitrary
polynomial function to the singular values of a matrix, known as Quantum
Singular Value Transformation (QSVT) [14]. QSVT is quite general: many quantum
algorithms can be recast to this framework, and for several problems, better
quantum algorithms can be obtained [14, 34]. In particular, QSVT has been
extremely useful in obtaining optimal quantum algorithms for linear algebra.
For instance, using QSVT, given the block-encoding of a matrix $A$, one could
obtain $A^{-c}$ with $c\in[0,\infty)$ with optimal complexity and by using
fewer additional qubits than prior art. This section briefly describes this
framework, which is a generalization of Quantum Signal Processing (QSP) [26,
Section 2], [25, Theorem 2], [32]. The reader may refer to [34] for a more
pedagogical overview of these techniques.

Let us begin by discussing the framework of Quantum Signal Processing. QSP is
a quantum algorithm to apply a $d$-degree bounded polynomial transformation
with parity $d\mod 2$ to an arbitrary quantum subsystem, using a quantum
circuit $U_{\Phi}$ consisting of only controlled single qubit rotations. This
is achieved by interleaving a signal rotation operator $W$ (which is an
$x$-rotation by some fixed angle $\theta$) and a signal processing operator
$S_{\phi}$ (which is a $z$-rotation by a variable angle $\phi\in[0,2\pi]$). In
this formulation, the signal rotation operator is defined as

$W(x):=\begin{pmatrix}x&i\sqrt{1-x^{2}}\\\ i\sqrt{1-x^{2}}&x\end{pmatrix},$
(12)

which is an $x$-rotation by angle $\theta=-2\arccos(x)$, and the signal
processing operator is defined as

$S_{\phi}:=e^{i\phi Z},$ (13)

which is a $z$-rotation by an angle $-2\phi$. Interestingly, sandwiching them
together for some
$\Phi:=(\phi_{0},\phi_{1},\ldots\phi_{d})\in\mathbb{R}^{d+1}$, as shown in
Equation 14, gives us a matrix whose elements are polynomial transformations
of $x$,

$\displaystyle U_{\Phi}$
$\displaystyle:=e^{i\phi_{0}Z}\prod_{j=1}^{j=d}\left(W(x)e^{i\phi_{j}Z}\right)$
(14) $\displaystyle=\begin{pmatrix}P(x)&iQ(x)\sqrt{1-x^{2}}\\\
iQ^{*}(x)\sqrt{1-x^{2}}&P^{*}(x)\end{pmatrix},$ (15)

such that

  1. 1.

$\deg P\leq d;\ \deg Q\leq d-1$,

  2. 2.

$P(x)$ has a parity $d\mod 2$,

  3. 3.

$|P(x)|^{2}+(1-x^{2})|Q(x)|^{2}=1\quad\forall x\in[-1,1]$.

Following the application of the quantum circuit $U_{\Phi}$ for an appropriate
$\Phi$, one can project into the top left block of $U_{\Phi}$ to recover the
polynomial $\bra{0}U_{\Phi}\ket{0}=P(x)$. Projecting to other basis allows the
ability to perform more interesting polynomial transformations, which can be
linear combinations of $P(x),Q(x)$, and their complex conjugates. For example,
projecting to $\\{\ket{+},\ket{-}\\}$ basis gives us

$\bra{+}U_{\Phi}\ket{+}=(P(x))+i(Q(x))\sqrt{1-x^{2}}.$ (16)

Quantum Signal Processing can be formally stated as follows.

######  Theorem 5 (Quantum Signal Processing, Corollary 8 from [14]).

Let $P\in\mathbb{C}[x]$ be a polynomial of degree $d\geq 2$, such that

  * •

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

  * •

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

  * •

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

  * •

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

Then there exists a $\Phi\in\mathbb{R}^{d}$ such that

$\prod_{j=1}^{d}\left(e^{i\phi_{j}\sigma_{z}}W(x)\right)=\begin{pmatrix}P(x)&\cdot\\\
\cdot&\cdot\end{pmatrix}.$ (17)

Thus, QSP allows us to implement any polynomial $P(x)$ that satisfies the
aforementioned requirements. Throughout this article, we refer to any such
polynomial $P(x)$ as a QSP polynomial. Quantum Singular Value Transformation
is a natural generalization of this procedure. It allows us to apply a QSP
polynomial transformation to each singular value of an arbitrary block of a
unitary matrix. In addition to this generalization, QSVT relies on the
observation that several functions can be well-approximated by QSP
polynomials. Thus, through QSVT one can transform each singular value of a
block-encoded matrix by any function that can be approximated by a QSP
polynomial. Since several linear algebra problems boil down to applying
specific transformations to the singular values of a matrix, QSVT is
particularly useful for developing fast algorithms for quantum linear algebra.
Next, we introduce QSVT formally via the following theorem.

######  Theorem 6 (Quantum Singular Value Transformation [13], Section 3.2).

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

$U_{\Phi}=\begin{cases}e^{i\phi_{1}(2\widetilde{\Pi}-I)}U_{A}\left[\prod_{k=1}^{(n-1)/2}e^{i\phi_{2k}(2\widetilde{\Pi}-I)}U_{A}^{\dagger}e^{i\phi_{2k+1}(2\widetilde{\Pi}-I)}U_{A}\right],&n\text{
is odd }\\\
\left[\prod_{k=1}^{n/2}e^{i\phi_{2k-1}(2\widetilde{\Pi}-I)}U_{A}^{\dagger}e^{i\phi_{2k}(2\widetilde{\Pi}-I)}U_{A}\right],&n\text{
is even},\end{cases}$ (18)

such that

$P^{SV}(A)=\begin{cases}\widetilde{\Pi}U_{\Phi}\Pi,&n\text{ is odd}\\\ \Pi
U_{\Phi}\Pi,&n\text{ is even},\end{cases}$ (19)

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

$P^{SV}(A):=\begin{cases}\sum_{j}P(\sigma_{j})\ket{v_{j}}\bra{w_{j}},&P\text{
is odd}\\\ \sum_{j}P(\sigma_{j})\ket{w_{j}}\bra{w_{j}},&P\text{ is
even}.\end{cases}$ (20)

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

$f(A):=\begin{cases}\sum_{j}f(\sigma_{j})\ket{v_{j}}\bra{w_{j}},&P\text{ is
odd}\\\ \sum_{j}f(\sigma_{j})\ket{w_{j}}\bra{w_{j}},&P\text{ is
even}.\end{cases}$ (21)

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

######  Lemma 7 (Robustness of Quantum Singular Value Transformation, [13],
Lemma 23).

Let $P\in\mathbb{C}[x]$ be a QSP polynomial of degree $n$. Let
$A,\tilde{A}\in\mathbb{C}^{N\times d}$ be matrices of spectral norm at most 1,
such that

$\norm{A-\tilde{A}}+\norm{\frac{A+\tilde{A}}{2}}^{2}\leq 1.$

Then,

$\norm{P^{SV}(A)-P^{SV}(\tilde{A})}\leq
n\sqrt{\frac{2}{1-\norm{\frac{A+\tilde{A}}{2}}^{2}}}\norm{A-\tilde{A}}.$

We will apply this theorem to develop a robust version of QSVT. More
precisely, in order to implement QSVT, we require access to a unitary $U_{A}$,
which is a block-encoding of some matrix $A$. This block-encoding, in most
practical scenarios, is not perfect: we only have access to a
$\varepsilon$-approximate block-encoding of $A$. If we want an
$\delta$-accurate implementation of $P^{SV}(A)$, how precise should the block-
encoding of $A$ be? Such a robustness analysis has been absent from prior work
involving QSVT and will allow us to develop robust versions of a number of
quantum algorithms in subsequent sections. The following theorem determines
the precision $\varepsilon$ required in the block-encoding of $A$ in terms of
$n$, the degree of the QSP polynomial that we wish to implement and $\delta$,
the accuracy of $P^{SV}(A)$.

######  Theorem 8 (Robust QSVT).

Let $P\in\mathbb{C}[x]$ be a QSP polynomial of degree $n\geq 2$. Let
$\delta\in[0,1]$ be the precision parameter. Let $U$ be an
$(\alpha,a,\varepsilon)$-block-encoding of matrix $A\in\mathbb{C}^{N\times d}$
satisfying $\norm{A}\leq\alpha/2$, implemented in cost $T$ for some
$\varepsilon\leq{\alpha\delta}/{2n}$. Then we can construct a
$(1,a+1,\delta)$-block-encoding of $P(A/\alpha)$ in cost $\order{nT}$.

###### Proof.

Let $\tilde{A}$ be the encoded block of $U$, then
$\norm{A-\tilde{A}}\leq\varepsilon$. Applying QSVT on $U$ with the polynomial
$P$, we get a block-encoding for $P(\tilde{A}/\alpha)$, with $\order{n}$ uses
of $U,U^{\dagger}$, and as many multiply-controlled NOT gates. Observe that
$\norm{\frac{A}{\alpha}-\frac{\tilde{A}}{\alpha}}\leq\frac{\varepsilon}{\alpha}\leq\frac{\delta}{2n}\leq\frac{1}{4}$,
and,

$\displaystyle\norm{\frac{\frac{A}{\alpha}+\frac{\tilde{A}}{\alpha}}{2}}^{2}=\norm{\frac{A}{\alpha}+\frac{\tilde{A}-A}{2\alpha}}^{2}\leq\left(\frac{\norm{A}}{\alpha}+\frac{\norm{\tilde{A}-A}}{2\alpha}\right)^{2}\leq\left(\frac{1}{2}+\frac{1}{8}\right)^{2}\leq\frac{1}{2}$

Therefore the error in the final block-encoding is given by invoking Lemma 7
with matrices $A/\alpha,\tilde{A}/\alpha$:

$\displaystyle\norm{P\left(\frac{A}{\alpha}\right)-P\left(\frac{\tilde{A}}{\alpha}\right)}\leq
n\sqrt{\frac{2}{1-\frac{1}{2}}}\leavevmode\nobreak\
\frac{\varepsilon}{\alpha}=\frac{2n\varepsilon}{\alpha}\leq\delta.$

∎

In Section 3, we will make use of Theorem 8, to develop robust quantum
algorithms for singular value discrimination, variable-time matrix inversion,
positive and negative powers of matrices. Subsequently, in Sec. 4, we shall
combine algorithmic primitives to design robust quantum regularized least
squares algorithms.

###  2.4 Variable Time Amplitude Amplification

Ambainis [3] defined the notion of a variable-stopping-time quantum algorithm
and formulated the technique of Variable Time Amplitude Amplification (VTAA),
a tool that can be used to amplify the success probability of a variable-
stopping-time quantum algorithm to a constant by taking advantage of the fact
that computation on some parts of an algorithm can complete earlier than on
other parts. The key idea here is to look at a quantum algorithm $\mathcal{A}$
acting on a state $\ket{\psi}$ as a combination of $m$ quantum sub-algorithms
$\mathcal{A}=\mathcal{A}_{m}\cdot\mathcal{A}_{m-1}\cdot\ldots\mathcal{A}_{1}$,
each acting on $\ket{\psi}$ conditioned on some ancilla flag being set.
Formally, a variable stopping time algorithm is defined as follows

######  Definition 9 (Variable-stopping-time Algorithm, [3]).

A quantum algorithm $\mathcal{A}$ acting on $\mathcal{H}$ that can be written
as $m$ quantum sub-algorithms,
$\mathcal{A}=\mathcal{A}_{m}\cdot\mathcal{A}_{m-1}\cdot\ldots\mathcal{A}_{1}$
is called a variable stopping time algorithm if
$\mathcal{H}=\mathcal{H}_{C}\otimes\mathcal{H}_{\mathcal{A}}$, where
$\mathcal{H}_{C}=\otimes_{i=1}^{m}\mathcal{H}_{C_{i}}$ with
$\mathcal{H}_{C_{i}}=\mathrm{span}(\ket{0},\ket{1})$, and each unitary
$\mathcal{A}_{j}$ acts on
$\mathcal{H}_{C_{j}}\otimes\mathcal{H}_{\mathcal{A}}$ controlled on the first
$j-1$ qubits $\ket{0}^{\otimes j-1}\in\otimes_{i=1}^{j-1}\mathcal{H}_{C_{i}}$
being in the all zero state.

Here $\mathcal{H}_{C_{i}}$ is a single qubit clock register. In VTAA,
$\mathcal{H}_{\mathcal{A}}$ has a flag space consisting of a single qubit to
indicate success,
$\mathcal{H}_{\mathcal{A}}=\mathcal{H}_{F}\otimes\mathcal{H}_{W}$. Here
$\mathcal{H}_{F}=\mathrm{Span}(\ket{g},\ket{b})$ flags the good and bad parts
of the run. Furthermore, for $1\leq i\leq m$, define the stopping times
$t_{i}$ such that $t_{1}<t_{2}<\cdots t_{m}=T_{\max}$, such that the algorithm
$\mathcal{A}_{j}\mathcal{A}_{j-1}\cdots\mathcal{A}_{1}$ having (gate/query)
complexity $t_{i}$ halts with probability

$p_{j}=\norm{\Pi_{C_{j}}\mathcal{A}_{j}\mathcal{A}_{j-1}\cdots\mathcal{A}_{1}\ket{0}_{\mathcal{H}}}^{2},$

where $\ket{0}_{\mathcal{H}}\in\mathcal{H}$ is the all zero quantum state and
$\Pi_{C_{j}}$ is the projector onto $\ket{1}$ in $\mathcal{H}_{C_{j}}$. From
this one can define the average stopping time of the algorithm $\mathcal{A}$
defined as

$\norm{T}_{2}=\sqrt{\sum_{j=1}^{m}p_{j}t^{2}_{j}}.$

For a variable stopping time algorithm if the average stopping time
$\norm{T}_{2}$ is less than the maximum stopping time $T_{\max}$, VTAA can
amplify the success probability $(p_{\mathrm{succ}})$ much faster than
standard amplitude amplification. In this framework, the success probability
of $\mathcal{A}$ is given by

$p_{\mathrm{succ}}=\norm{\Pi_{F}\mathcal{A}_{m}\mathcal{A}_{m-1}\cdots\mathcal{A}_{1}\ket{0}_{\mathcal{H}}}^{2}$

While standard amplitude amplification requires time scaling as
$\order{T_{\max}/\sqrt{p_{\mathrm{succ}}}}$, the complexity of VTAA is more
involved. Following [7], we define the complexity of VTAA as follows.

######  Lemma 10 (Efficient variable time amplitude amplification [6], Theorem
23).

Let $U$ be a state preparation unitary such that $U\ket{0}^{\otimes
k}=\sqrt{p_{\mathrm{prep}}}\ket{0}\ket{\psi_{0}}+\sqrt{1-p_{\mathrm{prep}}}\ket{1}\ket{\psi_{1}}$
that has a query complexity $T_{U}$. And let
$\mathcal{A}=\mathcal{A}_{m}\mathcal{A}_{m-1}\cdots\mathcal{A}_{1}$ be a
variable stopping time quantum algorithm that we want to apply to the state
$\ket{\psi_{0}}$, with the following known bounds: $p_{\mathrm{prep}}\geq
p^{\prime}_{\mathrm{prep}}$ and $p_{\mathrm{succ}}\geq
p^{\prime}_{\mathrm{succ}}$. Define $T^{\prime}_{\max}:=2T_{\max}/t_{1}$ and

$Q:=\left(T_{\max}+\frac{T_{U}+k}{\sqrt{p_{\mathrm{prep}}}}\right)\sqrt{\log\left(T^{\prime}_{\max}\right)}+\frac{\left(\norm{T}_{2}+\frac{T_{U}+k}{\sqrt{p_{\mathrm{prep}}}}\right)\log\left(T^{\prime}_{\max}\right)}{\sqrt{p_{\mathrm{succ}}}}.$

Then with success probability $\geq 1-\delta$, we can create a variable-
stopping time algorithm $\mathcal{A}^{\prime}$ that prepares the state
$a\ket{0}\mathcal{A}^{\prime}\ket{\psi_{0}}+\sqrt{1-a^{2}}\ket{1}\ket{\psi_{\textrm{garbage}}}$,
such that $a=\Theta(1)$ is a constant and $\mathcal{A}^{\prime}$ has the
complexity $\order{Q}$.

One cannot simply replace standard amplitude amplification with VTAA to boost
the success probability of a quantum algorithm. A crucial task would be to
recast the underlying algorithm in the VTAA framework. We will be applying
VTAA to the quantum algorithm for matrix inversion by QSVT. So, first of all,
in order to apply VTAA to the algorithm must be first recast into a variable-
time stopping algorithm so that VTAA can be applied.

Originally, Ambainis [3] used VTAA to improve the running time of the HHL
algorithm from $\order{\kappa^{2}\log N}$ to $\order{\kappa\log^{3}\kappa\log
N}$. Childs et al. [9] designed a quantum linear systems algorithm with a
polylogarithmic dependence on the accuracy. Additionally, they recast their
algorithm into a framework where VTAA could be applied to obtain a linear
dependence on $\kappa$. Later Chakraborty et al. [7] modified Ambainis’ VTAA
algorithm to perform variable time amplitude estimation.

In this work, to design quantum algorithms for $\ell_{2}$-regularized linear
regression, we use a quantum algorithm for matrix inversion by QSVT. We recast
this algorithm in the framework of VTAA to achieve nearly linear dependence in
$\kappa$ (the effective condition number of the matrix to be inverted). Using
QSVT instead of controlled Hamiltonian simulation improves the complexity of
the overall matrix inversion algorithm (QSVT and VTAA) by a log factor. It
also reduces the number of additional qubits substantially. Furthermore, we
replace a gapped quantum phase estimation procedure with a more efficient
quantum singular value discrimination algorithm using QSVT. This further
reduces the number of additional qubits by $O(\log^{2}(\kappa/\delta))$ than
in Refs. [9, 7], where $\kappa$ is the condition number of the underlying
matrix and $\delta$ is the desired accuracy. The details of the variable
stopping time quantum algorithm for matrix inversion by QSVT are laid out in
Section 3.3.

##  3 Algorithmic Primitives

This section introduces the building blocks of our quantum algorithms for
quantum linear regression with general $\ell_{2}$-regularization. As mentioned
previously, we work in the block-encoding framework. We develop robust quantum
algorithms for arithmetic operations, inversion, and positive and negative
powers of matrices using quantum singular value transformation, assuming we
have access to approximate block-encodings of these matrices. While some of
these results were previously derived assuming perfect block-encodings [14,
7], we calculate the precision required in the input block-encodings to output
a block-encoding or quantum state arbitrarily close to the target.

Given a $(\alpha,a,\varepsilon)$-block-encoding of a matrix $A$, we can
efficiently amplify the sub-normalization factor from $\alpha$ to a constant
and obtain an amplified block-encoding of $A$. For our quantum algorithms in
Sec. 4, we show working with pre-amplified block-encodings often yields better
complexities. We state the following lemma which was proven in Ref. [24]:

######  Lemma 11 (Uniform Block Amplification of Contractions, [24]).

Let $A\in\mathbb{R}^{N\times d}$ such that $\norm{A}\leq 1$ If $\alpha\geq 1$
and $U$ is a $(\alpha,a,\varepsilon)$-block-encoding of A that can be
implemented at a cost of $T_{U}$, then there is a
$(\sqrt{2},a+1,\varepsilon+\gamma)$-block-encoding of A that can be
implemented at a cost of $\order{\alpha T_{U}\log\left(1/\gamma\right)}$.

######  Corollary 12 (Uniform Block Amplification).

Let $A\in\mathbb{R}^{N\times d}$ and $\delta\in(0,1]$. Suppose $U$ is a
$(\alpha,a,\varepsilon)$-block-encoding of A, such that
$\varepsilon\leq\frac{\delta}{2}$, that can be implemented at a cost of
$T_{U}$. Then a $(\sqrt{2}\norm{A},a+1,\delta)$-block-encoding of A can be
implemented at a cost of $\order{\frac{\alpha
T_{U}}{\norm{A}}\log\left(\norm{A}/\delta\right)}$.

We now obtain the complexity of applying a block-encoded matrix to a quantum
state, which is a generalization of a lemma proven in Ref. [7].

######  Lemma 13 (Applying a Block-encoded Matrix on a Quantum State).

Let $A$ be an $s$-qubit operator such that its non-zero singular values lie in
$[\norm{A}/\kappa,\norm{A}]$. Also let $\delta\in(0,1)$, and $U_{A}$ be an
$(\alpha,a,\varepsilon)$-block-encoding of $A$, implementable in time $T_{A}$,
such that

$\varepsilon\leq\frac{\delta\norm{A}}{2\kappa}.$

Furthermore, suppose $\ket{b}$ be an $s$-qubit quantum state, prepared in time
$T_{b}$. Then we can prepare a state that is $\delta$-close to
$\frac{A\ket{b}}{\norm{A\ket{b}}}$ with success probability
$\Omega\left(1\right)$ at a cost of

$\order{\frac{\alpha\kappa}{\norm{A}}(T_{A}+T_{b})}$

######  Corollary 14 (Applying a pre-amplified Block-encoded Matrix on a
Quantum State).

Let $A$ be an $s$-qubit operator such that its non-zero singular values lie in
$[\norm{A}/\kappa,\norm{A}]$. Also let $\delta\in(0,1)$, and $U_{A}$ be an
$(\alpha,a,\varepsilon)$-block-encoding of $A$, implementable in time $T_{A}$,
such that

$\varepsilon\leq\frac{\delta\norm{A}}{4\kappa}.$

Furthermore, suppose $\ket{b}$ be an $s$-qubit quantum state that can be
prepared in time $T_{b}$. Then we can prepare a state that is $\delta$-close
to $\frac{A\ket{b}}{\norm{A\ket{b}}}$ with success probability
$\Omega\left(1\right)$ at a cost of

$\order{\frac{\alpha\kappa}{\norm{A}}\log\left(\frac{\kappa}{\delta}\right)T_{A}+\kappa
T_{b}}$

Now, it may happen that $U_{b}$ prepares a quantum state that is only
$\varepsilon$-close to the desired state $\ket{b}$. In such cases, we have the
following lemma

######  Lemma 15 (Robustness of state preparation).

Let $A$ be an $s$-qubit operator such that its non-zero singular values lie in
$\left[\norm{A}/\kappa,\norm{A}\right]$. Suppose $\ket{b^{\prime}}$ is a
quantum state that is $\varepsilon/2\kappa$-close to $\ket{b}$ and
$\ket{\psi}$ is a quantum state that is $\varepsilon/2$-close to
$A\ket{b^{\prime}}/\norm{A\ket{b^{\prime}}}$. Then we have that $\ket{\psi}$
is $\varepsilon$-close to $A\ket{b}/\norm{A\ket{b}}$.

The proofs for 12, Lemma 13, 14, and Lemma 15 can be found in Appendix A.

###  3.1 Arithmetic with Block-Encoded Matrices

The block-encoding framework embeds a matrix on the top left block of a larger
unitary $U$. It has been demonstrated that this framework allows us to obtain
sums, products, linear combinations of block-encoded matrices. This is
particularly useful for solving linear algebra problems in general. Here, we
state some of the arithmetic operations on block-encoded matrices that we
shall be using in order to design the quantum algorithms of Section 4 and
tailor existing results to our requirements.

First we prove a slightly more general form of linear combination of unitaries
in the block-encoding framework, presented in [14]. To do this we assume that
we are given optimal state preparation pairs, defined as follows.

######  Definition 16 (Optimal State Preparation Unitary).

Let $m\in\mathbb{Z}^{+}$, and $s=\lceil\log{m}\rceil$. Let
$\eta\in\mathbb{R}^{m}_{+}$. Then we call a $s$-qubit unitary $P$ a $\eta$
state-preparation unitary if

$P\ket{0}=\frac{1}{\sum_{j}\eta_{j}}\sum_{j}\sqrt{\eta_{j}}\ket{j}$

######  Lemma 17 (Linear Combination of Block Encoded Matrices).

For each $j\in\\{0,\ldots,m-1\\}$, let $A_{j}$ be an $s$-qubit operator, and
$y_{j}\in\mathbb{R}^{+}$. Let $U_{j}$ be a
$(\alpha_{j},a_{j},\varepsilon_{j})$-block-encoding of $A_{j}$, implemented in
time $T_{j}$. Define the matrix $A=\sum_{j}y_{j}A_{j}$, and the vector
$\eta\in\mathbb{R}^{m}$ s.t. $\eta_{j}=y_{j}\alpha_{j}$. Let $U_{\eta}$ be a
$\eta$ state-preparation unitary, implemented in time $T_{\eta}$. Then we can
implement a

$\left(\sum_{j}y_{j}\alpha_{j},\max_{j}(a_{j})+s,\sum_{j}y_{j}\varepsilon_{j}\right)$

block-encoding of $A$ at a cost of $\order{\sum_{j}T_{j}+T_{\eta}}$.

The proof is similar to the one in Ref. [14], with some improvements to the
bounds. The detailed proof can be found in Appendix A. We now specialize the
above lemma for the case where we need a linear combination of just two
unitaries. This is the case used in this work, and we obtain a better error
scaling for this by giving an explicit state preparation unitary.

######  Corollary 18 (Linear Combination of Two Block Encoded Matrices).

For $j\in\\{0,1\\}$, let $A_{j}$ be an $s$-qubit operator and
$y_{j}\in\mathbb{R}^{+}$. Let $U_{j}$ be a
$(\alpha_{j},a_{j},\varepsilon_{j})$-block-encoding of $A_{j}$, implemented in
time $T_{j}$. Then we can implement a
$(y_{0}\alpha_{0}+y_{1}\alpha_{1},1+\max(a_{0},a_{1}),y_{0}\varepsilon_{0}+y_{1}\varepsilon_{1})$
encoding of $y_{0}A_{0}+y_{1}A_{1}$ in time $\order{T_{0}+T_{1}}$.

###### Proof.

Let $\alpha=y_{0}\alpha_{0}+y_{1}\alpha_{1}$ and
$P=\frac{1}{\sqrt{\alpha}}\begin{pmatrix}\sqrt{y_{0}\alpha_{0}}&-\sqrt{y_{1}\alpha_{1}}\\\
\sqrt{y_{1}\alpha_{1}}&\sqrt{y_{0}\alpha_{0}}\end{pmatrix}$. By Definition 16,
$P$ is a $\\{y_{0}\alpha_{0},y_{1}\alpha_{1}\\}$ state-preparation-unitary.
Invoking Lemma 17 with $P$, we get the required unitary.

∎

Given block-encodings of two matrices $A$ and $B$, it is easy to obtain a
block-encoding of $AB$.

######  Lemma 19 (Product of Block Encodings, [13], Lemma 53).

If $U_{A}$ is an $(\alpha,a,\delta)$-block-encoding of an $s$-qubit operator
$A$ implemented in time $T_{A}$, and $U_{B}$ is a
$(\beta,b,\varepsilon)$-block-encoding of an $s$-qubit operator $B$
implemented in time $T_{B}$, then $(I^{\otimes b}\otimes U_{A})(I^{\otimes
a}\otimes U_{B})$ is an
$(\alpha\beta,a+b,\alpha\varepsilon+\beta\delta)$-block-encoding of $AB$
implemented at a cost of $\order{T_{A}+T_{B}}$.

Directly applying Lemma 19 results in a block-encoding of
$\frac{AB}{\alpha\beta}$. If $\alpha$ and $\beta$ are large, then the sub-
normalization factor $\alpha\beta$ might incur an undesirable overhead to the
cost of the algorithm that uses it. In many cases, the complexity of obtaining
products of block-encodings can be improved if we first amplify the block-
encodings (using Lemma 12) and then apply Lemma 19. We prove the following
lemma:

######  Lemma 20 (Product of Amplified Block-Encodings).

Let $\delta\in(0,1]$. If $U_{A}$ is an
$(\alpha_{A},a_{A},\varepsilon_{A})$-block-encoding of an $s$-qubit operator
$A$ implemented in time $T_{A}$, and $U_{B}$ is a
$(\alpha_{B},a_{B},\varepsilon_{B})$-block-encoding of an $s$-qubit operator
$B$ implemented in time $T_{B}$, such that
$\varepsilon_{A}\leq\frac{\delta}{4\sqrt{2}\norm{B}}$ and
$\varepsilon_{B}\leq\frac{\delta}{4\sqrt{2}\norm{A}}$. Then we can implement a
$(2\norm{A}\norm{B},a_{A}+a_{B}+2,\delta)$-block-encoding of $AB$ implemented
at a cost of

$\order{\left(\frac{\alpha_{A}}{\norm{A}}T_{A}+\frac{\alpha_{B}}{\norm{B}}T_{B}\right)\log\left(\frac{\norm{A}\norm{B}}{\delta}\right)}.$

###### Proof.

Using 12 for some $\delta_{A}\geq 2\varepsilon_{A}$ we get a
$(\sqrt{2}\norm{A},a_{A}+1,\delta_{A})$-block-encoding of $A$ at a cost of

$\order{\frac{\alpha_{A}T_{A}}{\norm{A}}\log\left(\norm{A}/\delta_{A}\right)}.$

Similarly for some $\delta_{B}\geq 2\varepsilon_{B}$ we get a
$(\sqrt{2}\norm{B},a_{B}+1,\delta_{B})$-block-encoding of $B$ at a cost of

$\order{\frac{\alpha_{B}T_{B}}{\norm{B}}\log\left(\norm{B}/\delta_{B}\right)}.$

Now using Lemma 19 we get a
$(2,a_{A}+a_{B}+2,\sqrt{2}\left(\norm{A}\delta_{B}+\norm{B}\delta_{A}\right))$-block-
encoding of $AB$. We can choose $\delta_{A}=\frac{\delta}{2\sqrt{2}\norm{B}}$
and $\delta_{B}=\frac{\delta}{2\sqrt{2}\norm{A}}$ which bounds the final
block-encoding error by $\delta$. ∎

Observe that we have assumed that $A$ and $B$ are $s$-qubit operators. For any
two matrices of dimension ${N\times d}$ and $d\times K$, such that $N,d,K\leq
2^{s}$, we can always pad them with rows and columns of zero entries and
convert them to $s$-qubit operators. Thus, in the scenario where $A$ and $B$
are not $s$-qubit operators, one can consider block encodings of padded
versions of these matrices. Note that this does not affect the operations on
the sub-matrix blocks encoding $A$ and $B$. Thus, the above results can be
used to perform block-encoded matrix products for arbitrary (compatible)
matrices.

Next we show how to find the block encoding of tensor product of matrices from
their block encodings. This procedure will be useful in creating the dilated
matrices required for regularization. The proof can be found in Appendix A.

######  Lemma 21 (Tensor Product of Block Encoded Matrices).

Let $U_{1}$ and $U_{2}$ be $(\alpha,a,\varepsilon_{1})$ and
$(\beta,b,\varepsilon_{2})$-block-encodings of $A_{1}$ and $A_{2}$, $s$ and
$t$-qubit operators, implemented in time $T_{1}$ and $T_{2}$ respectively.
Define $S:=\Pi_{i=1}^{s}\texttt{SWAP}_{a+b+i}^{a+i}$. Then, $S(U_{1}\otimes
U_{2})S^{\dagger}$ is an
$(\alpha\beta,a+b,\alpha\varepsilon_{2}+\beta\varepsilon_{1}+\varepsilon_{1}\varepsilon_{2})$
block-encoding of $A_{1}\otimes A_{2}$, implemented at a cost of
$\order{T_{1}+T_{2}}$.

We will now use Lemma 21 to augment one matrix into another, given their
approximate block-encodings.

######  Lemma 22 (Block-encoding of augmented matrix).

If $U_{A}$ is an $(\alpha_{A},a_{A},\varepsilon_{A})$-block encoding of an
$s$-qubit operator $A$ that can be implemented in time $T_{A}$ and $U_{B}$ is
an $(\alpha_{B},a_{B},\varepsilon_{B})$-block encoding of an $s$-qubit
operator $B$ that can be implemented in time $T_{B}$, then we an implement an
$(\alpha_{A}+\alpha_{B},\max(a_{A},a_{B})+2,\varepsilon_{A}+\varepsilon_{B})$-block-
encoding of

$A_{B}=\begin{pmatrix}A&0\\\ B&0\end{pmatrix}$

at a cost of $\order{T_{A}+T_{B}}$.

###### Proof.

Let $M_{A}=\begin{pmatrix}1&0\\\ 0&0\end{pmatrix}$. Then the SWAP gate is a
$(1,1,0)$ block encoding of $M_{A}$. By Lemma 21, we can implement
$U_{A}^{\prime}$, an $(\alpha_{A},a_{A}+1,\varepsilon_{A})$-block-encoding of
their tensor product $M_{A}\otimes A=\begin{pmatrix}A&0\\\ 0&0\end{pmatrix}$
at a cost of $\order{T_{A}}$. Similarly, Let $M_{B}=\begin{pmatrix}0&0\\\
1&0\end{pmatrix}$. Then $(I\otimes X)\cdot\texttt{SWAP}$ is a $(1,1,0)$-block-
encoding of $M_{B}$. Similarly Lemma 21, we can implement $U_{B}^{\prime}$, an
$(\alpha_{B},a_{B}+1,\varepsilon_{B})$-block-encoding of $M_{B}\otimes
B=\begin{pmatrix}0&0\\\ B&0\end{pmatrix}$ at a cost of $\order{T_{B}}$. We add
them by using 18 on $U_{A}^{\prime}$ and $U_{B}^{\prime}$, to implement
$U_{A_{B}}$, an
$(\alpha_{A}+\alpha_{B},2+\max(a_{A},a_{B}),\varepsilon_{A}+\varepsilon_{B})$-block-
encoding of $A_{B}=\begin{pmatrix}A&0\\\ B&0\end{pmatrix}$. This can be
implemented at a cost of $\order{T_{A}+T_{B}}$. ∎

###  3.2 Robust Quantum Singular Value Discrimination

The problem of deciding whether the eigenvalues of a Hamiltonian lie above or
below a certain threshold, known as eigenvalue discrimination, finds
widespread applications. For instance, the problem of determining whether the
ground energy of a generic local Hamiltonian is $<\lambda_{a}$ or
$>\lambda_{b}$ is known to be QMA-Complete [21]. Nevertheless, quantum
eigenvalue discrimination has been useful in preparing ground states of
Hamiltonians. Generally, a variant of quantum phase estimation, which
effectively performs a projection onto the eigenbasis of the underlying
Hamiltonian, is used to perform eigenvalue discrimination [16]. Recently, it
has been shown that QSVT can be used to approximate a projection onto the
eigenspace of an operator by implementing a polynomial approximation of the
sign function [29]. This was then used to design improved quantum algorithms
for ground state preparation.

In our work, we design a more general primitive, known as Quantum Singular
Value Discrimination (QSVD). Instead of eigenvalues, the algorithm
distinguishes whether a singular value $\sigma$ is $\leq\sigma_{a}$ or
$\geq\sigma_{b}$. This is particularly useful when the block-encoded matrix is
not necessarily Hermitian and hence, may not have well-defined eigenvalues. We
use this procedure to develop a more space-efficient variable stopping time
matrix inversion algorithm in Section 3.3. Owing to the widespread use of
singular values in a plethora of fields, we believe that our QSVD procedure is
of independent interest.

Let us define the sign function $\mathrm{sign}:\mathbb{R}\to\mathbb{R}$ as
follows:

$\mathrm{sign}(x)=\begin{cases}-1&\quad x<0\\\ 0&\quad x=0\\\ 1&\quad x>0.\\\
\end{cases}$ (22)

Given a threshold singular value $c$, Low and Chuang [24] showed that there
exists a polynomial approximation to $\mathrm{sign}(c-x)$ (based on its
approximation of the erf function). We use the result of Ref. [34], where such
a polynomial of even parity was considered. This is crucial, as for even
polynomials, QSVT maps right (left) singular vectors to right (left) singular
vectors, which enables us to use the polynomial in [34] for singular value
discrimination.

######  Lemma 23 (Polynomial approximation to the sign function [24, 28, 34]).

For any $\varepsilon,\Delta,c\in(0,1)$, there exists an efficiently computable
even polynomial $P_{\varepsilon,\Delta,c}(x)$ of degree
$l=\order{\frac{1}{\Delta}\log(1/\varepsilon)}$ such that

  * 1\.

$\forall x\in[0,1]\colon\absolutevalue{P_{\varepsilon,\Delta,c}(x)}\leq 1$

  * 2\.

$\forall
x\in[0,1]\setminus\left(c-\frac{\Delta}{2},c+\frac{\Delta}{2}\right)\colon\absolutevalue{P_{\varepsilon,\Delta,c}(x)-\mathrm{sign}(c-x)}\leq\varepsilon$

Therefore, given a matrix $A$ with singular values between $[0,1]$, we can use
QSVT to implement $P_{\varepsilon,\Delta,c}(A)$ which correctly distinguishes
between singular values of $A$ whose value is less than $c-\Delta/2$ and those
whose value is greater than $c+\Delta/2$. For our purposes, we shall consider
that we are given $U_{A}$, which is an $(\alpha,a,\varepsilon)$ block-encoding
of a matrix $A$. Our goal would be to distinguish whether a certain singular
value $\sigma$ satisfies $0\leq\sigma\leq\varphi$ or $2\varphi\leq\sigma\leq
1$. Since $U_{A}$ (approximately) implements $A/\alpha$, the task can be
rephrased as distinguishing whether a singular value of $A/\alpha$ is in
$[0,\varphi/\alpha]$ or in $[2\varphi/\alpha,1]$. For this, we develop a
robust version of quantum singular value discrimination ($QSVD(\phi,\delta)$),
which indicates the precision $\varepsilon$ required to commit an error that
is at most $\delta$.

######  Theorem 24 (Quantum Singular Value Discrimination using QSVT).

Suppose $A\in\mathbb{C}^{N\times N}$ is an $s$-qubit operator (where
$N=2^{s}$) with singular value decomposition
$A=\sum_{j\in[N]}\sigma_{j}\outerproduct{u_{j}}{v_{j}}$ such that all
$\sigma_{j}$ lie in the range $[0,1]$. Let
$\varphi\in\left(0,\frac{1}{2}\right)$ and $\delta\in(0,1]$ be some
parameters. Suppose that for some $\alpha\geq 2$ and $\varepsilon$ satisfying

$\varepsilon=o\left(\frac{\delta\varphi}{\log(1/\delta)}\right)$

we have access to $U_{A}$, an $(\alpha,a,\varepsilon)$-block-encoding of $A$
implemented in cost $T_{A}$. Then there exists a quantum algorithm
$QSVD(\varphi,\delta)$ which implements a $(1,a+1,\delta)$-block-encoding of
some $(s+1)$-qubit operator $D\in\mathbb{C}^{2N\times 2N}$ satisfying the
following constraints for all $j\in[N]$:

  * •

$\sigma_{j}\leq\varphi\implies D\ket{0}\ket{v_{j}}=\ket{0}\ket{v_{j}}$

  * •

$\sigma_{j}\geq 2\varphi\implies D\ket{0}\ket{v_{j}}=\ket{1}\ket{v_{j}}$

This algorithm has a cost of

$\order{\frac{\alpha}{\varphi}\log\left(\frac{1}{\delta}\right)T_{A}}.$

###### Proof.

We invoke Lemma 23 with parameters $\varepsilon^{\prime}:=\frac{\delta}{2}$,
$c:=\frac{3\varphi}{2\alpha}$ and $\Delta:=\frac{\varphi}{2\alpha}$, to
construct an even polynomial $P:=P_{\varepsilon^{\prime},\Delta,c}$ of degree
$n:=\order{\frac{\alpha}{\varphi}\log(\frac{1}{\varepsilon^{\prime}})}$, which
is an $\varepsilon^{\prime}$-approximation of
$f(x):=\text{sign}\left(\frac{3\varphi}{2\alpha}-x\right)$ for
$x\in\left[0,\frac{\varphi}{\alpha}\right]\cup\left[\frac{2\varphi}{\alpha},1\right]$.
Invoking Theorem 8 with $P$ and $U_{A}$, we get $U_{B}$ – a
$(1,a+1,\gamma)$-block-encoding of $B:=P(A/\alpha)$, implemented in cost
$\order{nT_{A}}$, where $\varepsilon$ must satisfy
$\varepsilon\leq\alpha\gamma/2n$.

Now consider the following unitary $W$ that acts on $s+a+2$ qubits:

$W:=\texttt{SWAP}_{[s,s+a+1]}^{\dagger}(H\otimes
I_{s+a+1})\left(C\text{-}U_{B}\right)(H\otimes
I_{s+a+1})\texttt{SWAP}_{[s,s+a+1]}$

$W$ is the required block-encoding of $D$, and $\texttt{SWAP}_{[l,r]}$
sequentially swaps adjacent qubits with indices in range $[l,r]$ effectively
moving qubit indexed $l$ to the right of qubit $r$. (where qubits are zero-
indexed, with higher indices for ancillas). Let $\tilde{B}$ be the top-left
block of $U_{B}$ (therefore $\norm{B-\tilde{B}}\leq\gamma$). Then we can
extract $\tilde{D}$, the top-left block of $W$ as follows:

$\displaystyle\tilde{D}$ $\displaystyle=\left(\bra{0}^{\otimes a+1}\otimes
I_{s+1}\right)\texttt{SWAP}_{[s,s+a+1]}^{\dagger}\left(\outerproduct{+}{+}\otimes
I_{s+a+1}+\outerproduct{-}{-}\otimes
U_{B}\right)\texttt{SWAP}_{[s,s+a+1]}\left(\ket{0}^{\otimes a+1}\otimes
I_{s+1}\right)$ $\displaystyle=\outerproduct{+}{+}\otimes
I_{s}+\outerproduct{-}{-}\otimes\tilde{B}$

Let us define index sets $L,R\subseteq[N]$ where
$L:=\\{j\in[N]\mid\sigma_{j}\leq\varphi\\}$ and
$R:=\\{j\in[N]\mid\sigma_{j}\geq 2\varphi\\}$; and the corresponding subspace
projections $\Pi_{L}:=\sum_{j\in L}\outerproduct{v_{j}}{v_{j}}$,
$\Pi_{R}:=\sum_{j\in R}\outerproduct{v_{j}}{v_{j}}$, and
$\Pi_{\perp}:=I_{s}-\Pi_{L}-\Pi_{R}$. Using these we pick our required
operator $D$ as follows:

$D:=I\otimes\Pi_{L}+X\otimes\Pi_{R}+\tilde{D}(I\otimes\Pi_{\perp})$

That is, $D$ behaves as expected on the required subspace, and acts identical
to $\tilde{D}$ on the remaining space. The error in the block-encoding can be
computed as

$\displaystyle\norm{D-\tilde{D}}$
$\displaystyle=\norm{I\otimes\Pi_{L}+X\otimes\Pi_{R}+\tilde{D}(I\otimes\Pi_{\perp})-\tilde{D}}$
$\displaystyle=\norm{I\otimes\Pi_{L}+X\otimes\Pi_{R}-\tilde{D}(I\otimes(\Pi_{L}+\Pi_{R}))}$
$\displaystyle=\norm{(I\otimes I_{s}-\tilde{D})(I\otimes\Pi_{L})+(X\otimes
I_{s}-\tilde{D})(I\otimes\Pi_{R})}$
$\displaystyle=\norm{\left(\outerproduct{-}{-}\otimes(I_{s}-\tilde{B})\right)(I\otimes\Pi_{L})-\left(\outerproduct{-}{-}\otimes(I_{s}+\tilde{B})\right)(I\otimes\Pi_{R})}$
$\displaystyle=\norm{(I_{s}-\tilde{B})\Pi_{L}-(I_{s}+\tilde{B})\Pi_{R}}$
$\displaystyle=\norm{(I_{s}-B)\Pi_{L}-(I_{s}+B)\Pi_{R}+(B-\tilde{B})(\Pi_{L}-\Pi_{R})}$
$\displaystyle\leq\norm{(I_{s}-P(A/\alpha))\Pi_{L}-(I_{s}+P(A/\alpha))\Pi_{R}}+\norm{B-\tilde{B}}\norm{\Pi_{L}-\Pi_{R}}$
$\displaystyle\leq\varepsilon^{\prime}+\gamma$

We can choose $\gamma=\delta/2$, therefore

$\varepsilon\leq\frac{\alpha\delta}{4n}=o\left(\frac{\delta\varphi}{\log(\frac{1}{\delta})}\right)$

∎

In Section 3.3, we develop a variable stopping time quantum algorithm for
matrix inversion using QSVT. In order to recast the usual matrix inversion to
the VTAA framework, we need to be able to apply this algorithm to specific
ranges of the singular values of the matrix. This is achieved by applying a
controlled QSVD algorithm, to determine whether the input singular vector
corresponds to an singular value less than (or greater than) a certain
threshold. Based on the outcome of controlled QSVD, the standard inversion
algorithm is applied. These two steps correspond to sub-algorithms
$\mathcal{A}_{j}$ of the VTAA framework.

In prior works such as Refs. [3, 9, 7], gapped phase estimation (GPE) was used
to implement this. GPE requires an additional register of
$\order{\log(\kappa)\log(1/\delta)}$ qubits to store the estimated phases. For
the whole VTAA procedure, $\log\kappa$ such registers are needed. As a result,
substituting GPE with QSVD, we save $\order{\log^{2}(\kappa)\log(1/\delta)}$
qubits.

###  3.3 Variable-Time Quantum Algorithm for Matrix Inversion using QSVT

Matrix inversion by QSVT applies a polynomial approximation of $f(x)=1/x$,
satisfying the constraints laid out in Section 2.3. Here, we make use of the
result of [34] to implement $A^{+}$. We adapt their result to the scenario
where we have an approximate block-encoding of $A$ as input. Finally, we
convert this to a variable stopping time quantum algorithm and apply VTAA to
obtain a linear dependence on the condition number of $A$.

######  Lemma 25 (Matrix Inversion polynomial (Appendix C of [34])).

Given $\kappa\geq 1,\varepsilon\in\mathbb{R}^{+}$, there exists an odd QSP
polynomial $P_{\varepsilon,\kappa}^{\text{MI}}$ of degree
$\order{\kappa\log(\kappa/\varepsilon)}$, which is an
$\frac{\varepsilon}{2\kappa}$ approximation of the function
$f(x)=\frac{1}{2\kappa x}$ in the range
$\mathcal{D}:=[-1,-\frac{1}{\kappa}]\cup[\frac{1}{\kappa},1]$. Also in this
range $P_{\varepsilon,\kappa}^{\text{MI}}$ is bounded from above by $1$, i.e.
$\forall
x\in\mathcal{D}:\left|P_{\varepsilon,\kappa}^{\text{MI}}(x)\right|\leq 1$.

######  Theorem 26 (Inverting Normalized Matrices using QSVT).

Let $A$ be a normalized matrix with non-zero singular values in the range
$[1/\kappa_{A},1]$ for some $\kappa_{A}\geq 1$. Let $\delta\in(0,1].$ For some
$\varepsilon=o\left(\frac{\delta}{\kappa_{A}^{2}\log\left({\kappa_{A}}/{\delta}\right)}\right)$
and $\alpha\geq 2$, let $U_{A}$ be an $(\alpha,a,\varepsilon)$-block-encoding
of $A$, implemented in time $T_{A}$. Then we can implement a
$(2\kappa_{A},a+1,\delta)$-block-encoding of $A^{+}$ at a cost of

$\order{\kappa_{A}\alpha\log\left(\frac{\kappa_{A}}{\delta}\right)T_{A}}.$

###### Proof.

We use the matrix inversion polynomial defined in Lemma 25,
$P:=P^{MI}_{\phi,\kappa}$ for this task, with $\kappa=\kappa_{A}\alpha$ and an
appropriate $\phi$. This has a degree of
$n:=\order{\kappa_{A}\alpha\log\left(\kappa_{A}\alpha/\phi\right)}$. We invoke
Theorem 8 to apply QSVT using the polynomial $P$ above, block-encoding
$U_{A}$, and an appropriate error parameter $\gamma$ such that
$\varepsilon\leq\alpha\gamma/2n$, to get the unitary $U$, a
$(1,a+1,\gamma)$-block-encoding of $P(A/\alpha)$. As $P$ is a
$(\phi/2\kappa)$-approximation of $f(x):=1/2\kappa x$, we have

$\displaystyle\norm{f(A/\alpha)-P(A/\alpha)}\leq\frac{\phi}{2\kappa},$

which implies $U$ is a $(1,a+1,\gamma+\phi/2\kappa)$-block-encoding of
$f(A/\alpha)$. And because $f(A/\alpha)=\frac{\alpha
A^{+}}{2\kappa}=A^{+}/2\kappa_{A}$, we can re-interpret $U$ as a
$(2\kappa_{A},a+1,2\kappa_{A}\gamma+\phi/\alpha)$-block-encoding of $A^{+}$.
Choosing $2\kappa_{A}\gamma=\phi/\alpha=\delta/2$, the final block-encoding
has an error of $\delta$. This gives us $\phi=\alpha\delta/2$ and
$\gamma=\delta/4\kappa_{A}$, and

$\varepsilon\leq\frac{\alpha\gamma}{2n}=\frac{\alpha\delta}{8\kappa_{A}n}=\order{\frac{\delta}{\kappa_{A}^{2}\log(\kappa_{A}/\delta)}}$

∎

Next, we design a map $W(\gamma,\delta)$ that uses QSVT to invert the singular
values of a matrix if they belong to a particular domain. This helps us recast
the usual matrix inversion algorithm as a variable-stopping-time algorithm and
will be a key subroutine for boosting the success probability of this
algorithm using VTAA. This procedure was also used in Refs. [9, 7] for the
quantum linear systems algorithms.

######  Theorem 27 (Efficient inversion of block-encoded matrix).

Let $A$ be a normalized matrix with non-zero singular values in the range
$\left[1/\kappa,1\right]$, for some $\kappa\geq 1$. Let
$\delta\in(0,1];\leavevmode\nobreak\ 0<\gamma\leq 1$. Let $U_{A}$ be an
$(\alpha,a,\varepsilon)$-block-encoding of $A$ implemented in time $T_{A}$,
such that $\alpha\geq 2$ and
$\varepsilon=o\left(\frac{\delta\gamma^{2}}{\log\left(\frac{1}{\delta\gamma}\right)}\right)$
. Then for any quantum state $\ket{b}$ that is spanned by the left singular
vectors of $A$ corresponding to the singular values in the range
$\left[\gamma,1\right]$, there exists a unitary $W(\gamma,\delta)$ that
implements

$W(\gamma,\delta):\ket{0}_{F}\ket{0}_{Q}\ket{b}_{I}\mapsto\dfrac{1}{a_{\max}}\ket{1}_{F}\ket{0}_{Q}f(A)\ket{b}_{I}+\ket{0}_{F}\ket{\perp}_{QI}$
(23)

where $a_{\max}=\order{\kappa_{A}}$ is a constant independent of $\gamma$,
$\ket{\perp}_{QI}$ is an unnormalized quantum state orthogonal to
$\ket{0}_{Q}$ and $\norm{f(A)\ket{b}-A^{+}\ket{b}}\leq\delta$. Here $F$ is a
1-qubit flag register, $Q$ is an $\alpha$-qubit ancilla register, and $I$ is
the $\log N$-qubit input register. This unitary has a cost

$\order{\frac{\alpha}{\gamma}\log\left(\frac{1}{\gamma\delta}\right)T_{A}}$
(24)

###### Proof.

Since we only need to invert the singular values in a particular range, we can
use the procedure in Theorem 26 with $\kappa_{A}$ modified to the restricted
range. That gives us the description of a quantum circuit
$\widetilde{W}(\gamma,\delta)$ that can implement the following map

$\widetilde{W}(\gamma,\delta):\ket{b}_{I}\ket{0}_{Q}\mapsto\frac{\gamma}{2}f(A)\ket{b}_{I}\ket{0}_{Q}+\ket{\perp}_{QI},$

where $\ket{\perp}$ is an unnomalized state with no component along
$\ket{0}_{Q}$. This has the same cost as Equation 24. Here
$\norm{f(A)\ket{\psi}-A^{+}\ket{\psi}}\leq\delta$ whenever $\ket{\psi}$ is a
unit vector in the span of the singular vectors of $A$ corresponding to the
singular values in $[\gamma,1]$. This follows from the sub-multiplicativity
property of the matrix-vector product.

Next, we must transform the amplitude of the good part of the state to
$\Theta({\kappa})$, independent of $\gamma$. To achieve this, we will have to
flag it with an ancillary qubit to use a controlled rotation to modify the
amplitude. Thus we add a single qubit $\ket{0}_{F}$ register and flip this
register controlled on register $Q$ being in the state $\ket{0}$ (the good
part). This gives us the transformation

$\widetilde{W}^{\prime}(\gamma,\delta):\ket{0}_{F}\ket{b}_{I}\ket{0}_{Q}\mapsto\frac{\gamma}{2}\ket{1}_{F}f(A)\ket{b}_{I}\ket{0}_{Q}+\ket{0}_{F}\ket{\perp}_{QI}$

Then we use a controlled rotation to replace the amplitude $\gamma/2$ with
some constant $a_{\max}^{-1}$ which is independent of $\gamma$, which is
achieved by introducing the relevant phase to the flag space

$\ket{1}_{F}\mapsto\frac{2}{\gamma
a_{\max}}\ket{1}_{F}+\sqrt{1-\frac{4}{\gamma^{2}a^{2}_{\max}}}\ket{0}_{F}.$

This gives us the desired $W(\gamma,\delta)$ as in Equation 23. ∎

Given such a unitary $W(\gamma,\delta)$, Ref. [7] laid out a procedure for a
variable time quantum algorithm $\mathcal{A}$ that takes as input the block
encoding of an $N\times d$ matrix $A$, and a state preparation procedure
$U_{b}:\ket{0}^{\otimes n}\mapsto\ket{b}$, and outputs a quantum state that is
a bounded distance away from $A^{+}\ket{b}/\norm{A^{+}\ket{b}}$. In order to
determine the branches of the algorithm on which to apply VTAA at a particular
iteration, [9, 7, 3] use the technique of gapped phase estimation, which given
a unitary $U$, a threshold $\phi$ and one of its eigenstate $\ket{\lambda}$,
decides if the corresponding eigenvalue is a bounded distance below the
threshold, or a bounded distance above it. In this work, we replace gapped
phase estimation with the QSVD algorithm (Theorem 24) which can be applied
directly to any block-encoded (not necessarily Hermitian) matrix $A$, and
allows for saving on $\order{\log^{2}(\kappa/\delta)}$ qubits.

The Variable time Algorithm: This algorithm will be a sequence of $m$ sub-
algorithms
$\mathcal{A}=\mathcal{A}_{m}\cdot\mathcal{A}_{m-1}\cdot\ldots\mathcal{A}_{1}$,
where $m=\lceil\log\kappa\rceil+1$. The overall algorithm acts on the
following registers:

  * •

$m$ single qubit clock registers $C_{i}:i\in[m]$.

  * •

An input register $I$, initialized to $\ket{0}^{\otimes s}$.

  * •

Ancillary register space $Q$ for the block encoding of $A$, initialized to
$\ket{0}^{\otimes a}$.

  * •

A single qubit flag register $\ket{0}_{F}$ used to flag success of the
algorithm.

Once we have prepared the above state space, we use the state preparation
procedure to prepare the state $\ket{b}$. Now we can define how each
$\mathcal{A}_{j}$ acts on the state space. Let
$\varepsilon^{\prime}=\frac{\delta}{a_{\max}m}$. The action of
$\mathcal{A}_{j}$ can be broken down into two parts:

  1. 1.

If $C_{j-1}\ldots C_{1}$ is in state $\ket{0}^{\otimes(j-1)}$, apply
QSVD$(2^{-j},\varepsilon^{\prime})$, (Theorem 24) to the state $\ket{b}$. The
output is to be written to the clock register $C_{j}$.

  2. 2.

If the state of $C_{j}$ is now $\ket{1}$, apply
$W(2^{-j},\varepsilon^{\prime})$ to $I\otimes F\otimes Q$.

Additionally, we would need algorithms
$\mathcal{A}^{\prime}=\mathcal{A}^{\prime}_{m}\cdots\mathcal{A}^{\prime}_{1}$
which are similar to $\mathcal{A}$, except that in Step 2, it implements
$W^{\prime}$ which sets the flag register to $1$. That is,

$W^{\prime}\ket{b}_{I}\ket{0}_{F}\ket{0}_{Q}=\ket{b}_{I}\ket{1}_{F}\ket{0}_{Q}.$

Now we are in a position to define the variable time quantum linear systems
algorithm using QSVT.

######  Theorem 28 (Variable Time Quantum Linear Systems Algorithm Using
QSVT).

Let $\varepsilon,\delta>0$. Let $A$ is a normalized $N\times d$ matrix such
that its non-zero singular values lie in $[1/\kappa,1]$. Suppose that for

$\varepsilon=o\left(\frac{\delta}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right),$

we have access to $U_{A}$ which is an $(\alpha,a,\varepsilon)$-block-encoding
of $A$, implemented with cost $T_{A}$. Let $\ket{b}$ be a state vector which
is spanned by the left singular vectors of $A$. Suppose there exists a
procedure to prepare the state $\ket{b}$ in cost $T_{b}$. Then there exists a
variable time quantum algorithm that outputs a state that is $\delta$-close
$\frac{A^{+}\ket{b}}{\norm{A^{+}\ket{b}}}$ at a cost of

$\order{\kappa\log\kappa\left(\alpha
T_{A}\log\left(\frac{\kappa}{\delta}\right)+T_{b}\right)}$ (25)

using $\order{\log\left(\kappa\right)}$ additional qubits.

###### Proof.

The correctness of the algorithm is similar to that of Refs. [9, 7], except
here, we use QSVD instead of gapped phase estimation. According to Lemma 10,
we need $T_{\max}$ (the maximum time any of the sub-algorithms
$\mathcal{A}_{j}$ take), $\norm{T}_{2}^{2}$ (the $\ell_{2}$-averaged stopping
time of the sub-algorithms), and $\sqrt{p_{\mathrm{succ}}}$ (the square root
of the success probability.) Now each sub-algorithm consists of two steps,
implementing QEVD with precision $2^{-j}$ and error $\varepsilon^{\prime}$,
followed by $W(2^{-j},\varepsilon^{\prime})$. From Theorem 24, the first step
costs

$\order{\alpha T_{A}2^{j}\log\left(\frac{1}{\varepsilon^{\prime}}\right)},$

and the cost of implementing $W(2^{-j},\varepsilon^{\prime})$ is as described
in Equation 24. Thus the overall cost of $\mathcal{A}_{j}$, which is the sum
of these two costs, turns out to be

$\order{\alpha T_{A}2^{j}\log\left(\frac{2^{j}}{\varepsilon^{\prime}}\right)}$
(26)

Note that the time $t_{j}$ required to implement
$\mathcal{A}_{j}\ldots\mathcal{A}_{1}$ is also the same as Equation 26. Also,

$\displaystyle T_{\max}$ $\displaystyle=\max_{j}\leavevmode\nobreak\ t_{j}$
$\displaystyle=\max_{j}\leavevmode\nobreak\ \order{\alpha
T_{A}2^{j}\log\left(\frac{2^{j}}{\varepsilon^{\prime}}\right)}$
$\displaystyle=\order{\alpha
T_{A}\kappa\log\left(\frac{\kappa}{\varepsilon^{\prime}}\right)}$
$\displaystyle=\order{\alpha
T_{A}\kappa\log\left(\frac{\kappa\log\left(\kappa\right)}{\delta}\right)}.$

The $\norm{T}_{2}^{2}$ is dependent on the probability that $\mathcal{A}$
stops at the $j^{\mathrm{th}}$ step. This is given by
$p_{j}=\norm{\Pi_{C_{j}}\mathcal{A}_{j}\ldots\mathcal{A}_{1}\ket{\psi}_{I}\ket{0}_{CFPQ}}^{2}$,
where $\Pi_{C_{j}}$ is the projector on $\ket{1}_{C_{j}}$, the
$j^{\mathrm{th}}$ clock register. From this, $\norm{T}_{2}^{2}$ can be
calculated as

$\displaystyle\norm{T}_{2}^{2}$ $\displaystyle=\sum_{j}p_{j}t_{j}^{2}$
$\displaystyle=\sum_{j}\norm{\Pi_{C_{j}}\mathcal{A}_{j}\ldots\mathcal{A}_{1}\ket{\psi}_{I}\ket{0}_{CFPQ}}^{2}t_{j}^{2}$
$\displaystyle=\sum_{k}\absolutevalue{c_{k}}^{2}\sum_{j}\left(\norm{\Pi_{C_{j}}\mathcal{A}_{j}\ldots\mathcal{A}_{1}\ket{v_{k}}_{I}\ket{0}_{CFPQ}}^{2}t_{j}^{2}\right)$
$\displaystyle=\order{\alpha^{2}T_{A}^{2}\sum_{k}{\log^{2}\left(\frac{1}{\sigma_{k}\varepsilon^{\prime}}\right)\frac{\absolutevalue{c_{k}}^{2}}{\sigma_{k}^{2}}}}$

Therefore

$\norm{T}_{2}=\order{\alpha
T_{A}\log\left(\frac{\kappa\log{\kappa}}{\delta}\right)\sqrt{\sum_{k}\frac{\absolutevalue{c_{k}}^{2}}{\sigma_{k}^{2}}}}.$
(27)

Next we calculate the success probability.

$\displaystyle\sqrt{p_{\mathrm{succ}}}$
$\displaystyle=\norm{\Pi_{F}\frac{A^{-1}}{\alpha_{\max}}\ket{b}_{I}\ket{\phi}_{CFPQ}}+\order{m\varepsilon^{\prime}}$
$\displaystyle=\frac{1}{\alpha_{\max}}\sqrt{\sum_{j}\frac{\absolutevalue{c_{j}}^{2}}{\sigma_{j}^{2}}}+\order{\frac{\delta}{\alpha_{\max}}}$
$\displaystyle=\Omega\left(\frac{1}{\kappa}\sqrt{\sum_{j}\frac{\absolutevalue{c_{j}}^{2}}{\sigma_{j}^{2}}}\right)$

Given these, we can use Lemma 10 to write the final complexity of matrix
inversion with VTAA:

$\displaystyle
T_{\max}+T_{b}+\frac{(\norm{T}_{2}+T_{b})\log\left(T^{\prime}_{\max}\right)}{\sqrt{p_{\textrm{succ}}}}=\order{\kappa\log\kappa\left(\alpha
T_{A}\log\left(\frac{\kappa}{\delta}\right)+T_{b}\right)}$

The upper bound on the precision required for the input block-encoding,
$\varepsilon$, can be calculated from the bounds on the precisions for
$W(\kappa,\varepsilon^{\prime})$ (Theorem 27) and
$\mathrm{QSVD}(\kappa,\varepsilon^{\prime})$ (Theorem 24) as follows:

$\displaystyle\varepsilon=o\left(\min\left(\frac{\varepsilon^{\prime}}{\kappa^{2}\log\left(\frac{\kappa}{\varepsilon^{\prime}}\right)},\frac{\varepsilon^{\prime}}{\kappa\log\left(\frac{1}{\varepsilon^{\prime}}\right)}\right)\right)=o\left({\frac{\varepsilon^{\prime}}{\kappa^{2}\log\left(\frac{\kappa}{\varepsilon^{\prime}}\right)}}\right)=o\left(\frac{\delta}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right)$

∎

The overall complexity is better by a $\log$ factor and requires
$\order{\log^{2}(\kappa/\delta)}$ fewer additional qubits as compared to the
variable time algorithms in Refs. [9, 7].

###  3.4 Negative Powers of Matrices using QSVT

We consider the problem: given an approximate block-encoding of a matrix $A$,
we need to prepare a block-encoding of $A^{-c}$, where $c\in(0,1)$. This
procedure will be used to develop algorithms for $\ell_{2}$-regularized
versions of GLS. We will directly use the results of [14].

######  Lemma 29 (Polynomial approximations of negative power functions, [13],
Corollary 67).

Let $\varepsilon,\delta\in(0,\frac{1}{2}],c>0$ and let
$f(x):=\frac{\delta^{c}}{2}x^{-c}$, then there exist even/odd polynomials
$P_{c,\varepsilon,\delta},P^{\prime}_{c,\varepsilon,\delta}\in\mathbb{R}[x]$
such that $\norm{P_{c,\varepsilon,\delta}-f}_{[\delta,1]}\leq\varepsilon$,
$\norm{P_{c,\varepsilon,\delta}}_{[-1,1]}\leq 1$ and
$\norm{P^{\prime}_{c,\varepsilon,\delta}-f}_{[\delta,1]}\leq\varepsilon$,
$\norm{P^{\prime}_{c,\varepsilon,\delta}}_{[-1,1]}\leq 1$. Moreover the degree
of the polynomials are
$\order{\frac{\max(1,c)}{\delta}\log\left(\frac{1}{\varepsilon}\right)}$.

######  Theorem 30 (Negative fractional powers of a normalized matrix using
QSVT).

Let $c\in(0,1)$ be some constant and $\delta\in(0,1]$ Let $A$ be a normalized
matrix with non-zero singular values in the range $[1/\kappa,1]$. Let $U_{A}$
be a $(\alpha,a,\varepsilon)$-block-encoding of a matrix $A$, implemented in
time $T_{A}$ such that $\alpha\geq 2$ and

$\varepsilon=o\left(\frac{\delta}{\kappa^{c+1}\log(\kappa/\delta)}\right)$

Then we can construct a $(2\kappa^{c},a+1,\delta)$-block-encoding of $A^{-c}$
at a cost of

$\order{\alpha\kappa\log\left(\frac{\kappa}{\delta}\right)\leavevmode\nobreak\
T_{A}}.$

###### Proof.

From Lemma 29, using $\Delta:=\frac{1}{\kappa\alpha}$ and an appropriate
$\varphi\in(0,\frac{1}{2}]$, we get an even QSP polynomial
$P:=P_{c,\varphi,\Delta}$ which is $\varphi$-close to
$f(x):=\frac{1}{2\kappa^{c}\alpha^{c}x^{c}}$, and has degree $n$ such that
$n=\order{\alpha\kappa\log\left(\frac{1}{\varphi}\right)}$. Therefore

$\norm{f(A/\alpha)-P(A/\alpha)}\leq\varphi.$

Using Theorem 8 we can construct $U_{P}$, a $(1,a+1,\gamma)$-block-encoding of
$P(A/\alpha)$, given that $\varepsilon\leq\frac{\alpha\gamma}{2n}$. Then from
triangle inequality it follows that it is a $(1,a+1,\leavevmode\nobreak\
\varphi+\gamma)$-block-encoding of $f(A/\alpha)$. And because
$f(A/\alpha)=\frac{A^{-c}}{2\kappa^{c}}$, $U_{P}$ can be re-interpreted as a
$(2\kappa^{c},a+1,2\kappa^{c}(\varphi+\gamma))$-block-encoding of $A^{-c}$. We
therefore choose $\varphi=\gamma=\frac{\delta}{4\kappa^{c}}$, and choose
$\varepsilon$ as

$\varepsilon=o\left(\alpha\frac{\delta}{4\kappa^{c}}\frac{1}{\alpha\kappa\log(4\kappa^{c}/\delta)}\right)=o\left(\frac{\delta}{\kappa^{c+1}\log(\kappa/\delta)}\right)$

∎

Having discussed the necessary algorithmic primitives, we are now in a
position to design quantum algorithms for linear regression with general
$\ell_{2}$-regularization. We will first deal with ordinary least squares
followed by weighted and generalized least squares.

##  4 Quantum Least Squares with General $\ell_{2}$-Regularization

In this section, we derive the main results of our paper, namely quantum
algorithms for quantum ordinary least squares (OLS), quantum weighted least
squares (WLS) and quantum generalized least squares (GLS) with
$\ell_{2}$-regularization.

###  4.1 Quantum Ordinary Least Squares

Given $N$ data points $\\{a_{i},b_{i}\\}_{i=1}^{N}$ such that
$a_{i}\in\mathbb{R}^{d}$ and $b_{i}\in\mathbb{R}$, the objective of linear
regression is to find $x\in\mathbb{R}^{d}$ that minimizes the loss function

$\mathcal{L}_{O}=\sum_{j=1}^{N}(x^{T}a_{i}-b_{i})^{2}.$ (28)

Consider the $N\times d$ matrix $A$ (known as the data matrix) such that the
$i^{\mathrm{th}}$ row of $A$ is the vector $a_{i}$ transposed and the column
vector $b=(b_{1}\cdots b_{N})^{T}$. Then, the solution to the OLS problem is
given by $x=(A^{T}A)^{-1}A^{T}b=A^{+}b$.

For the $\ell_{2}$-regularized version of the OLS problem, a penalty term is
added to its objective function. This has the effect of shrinking the singular
values of $A$ which helps overcome problems such as rank deficiency and
overfitting for the OLS problem. The loss function to be minimized is of the
form

$\norm{Ax-b}^{2}_{2}+\norm{Lx}^{2}_{2},$ (29)

where $L$ is the $N\times d$ penalty matrix and $\lambda>0$ is the optimal
regularizing parameter. The solution $x\in\mathbb{R}^{d}$ satisfies

$x=(A^{T}A+\lambda L^{T}L)^{-1}A^{T}b.$ (30)

Therefore, for quantum ordinary least squares with general
$\ell_{2}$-regularization, we assume that we have access to approximate block-
encodings of the data matrix $A$, $L$ and a procedure to prepare the quantum
state $\ket{b}=\sum_{j=1}^{N}b_{j}\ket{j}/\norm{b}$. Our algorithm outputs a
quantum state that is close to

$\ket{x}=\dfrac{(A^{T}A+\lambda
L^{T}L)^{-1}A^{T}\ket{b}}{\norm{(A^{T}A+\lambda L^{T}L)^{-1}A^{T}\ket{b}}}.$
(31)

In order to implement a quantum algorithm that implements this, a
straightforward approach would be the following: We first construct block-
encodings of $A^{T}A$ and $L^{T}L$, given block encodings of $A$ and $L$,
respectively (Using Lemma 19). We could then implement a block-encoding of
$A^{T}A+\lambda L^{T}L$ using these block encodings (By Lemma 17). On the
other hand, we could also prepare a quantum state proportional to
$A^{T}\ket{b}$ by using the block-encoding for $A$ and the unitary preparing
$\ket{b}$. Finally, using the block encoding of $A^{T}A+\lambda L^{T}L$, we
could implement a block-encoding of $(A^{T}A+\lambda L^{T}L)^{-1}$ (using
Theorem 26) and apply it to the state $A^{T}\ket{b}$. Although this procedure
would output a quantum state close to $\ket{x}$, it is not efficient. It is
easy to see that the inverse of $A^{T}A+\lambda L^{T}L$, would be implemented
with a complexity that has a quadratic dependence on the condition numbers of
$A$ and $L$. This would be undesirable as it would perform worse than the
unregularized quantum least squares algorithm, where one is able to implement
$A^{+}$ directly. However, it is possible to design a quantum algorithm that
performs significantly better than this.

The first observation is that it is possible to recast this problem as finding
the pseudoinverse of some augmented matrix. Given the data matrix
$A\in\mathbb{R}^{N\times d}$, the regularizing matrix $L\in\mathbb{R}^{N\times
d}$, let us define the following augmented matrix

$A_{L}:=\begin{pmatrix}A&0\\\ \sqrt{\lambda}L&0\end{pmatrix}.$ (32)

It is easy to see that the top left block of $A^{+}_{L}=(A^{T}A+\lambda
L^{T}L)^{-1}A^{T}$, which is the required linear transformation to be applied
to $b$. Consequently, our strategy would be to implement a block-encoding of
$A_{L}$, given block-encodings of $A$ and $L$. Following this, we use matrix
inversion by QSVT to implement $A^{+}_{L}\ket{b}\ket{0}$. The first register
is left in the quantum state given in Equation 31.

From this, it is clear that the complexity of our quantum algorithm would
depend on the effective condition number of the augmented matrix $A_{L}$. In
this regard, we shall assume that the penalty matrix $L$ is a good
regularizer. That is, $L$ is chosen such that it does not have zero singular
values (positive definite). This is a fair assumption as if $L$ has only non-
zero singular values, the minimum singular value of $A_{L}$ is guaranteed to
be lower bounded by the minimum singular value of $L$. This ensures that the
effective condition number of $A_{L}$ depends on $\kappa_{L}$, even when the
data matrix $A$ has zero singular values and $A^{T}A$ is not invertible.
Consequently, this also guarantees that regularized least squares provide an
advantage over their unregularized counterparts.

Next, we obtain bounds on the effective condition number of the augmented
matrix $A_{L}$ for a good regularizer $L$ via the following lemma:

######  Lemma 31 (Condition number and Spectral Norm of $A_{L}$).

Let the data matrix $A$ and the positive definite penalty matrix $L$ have
spectral norms $\norm{A}$ and $\norm{L}$, respectively. Furthermore, suppose
their effective condition numbers be upper bounded by $\kappa_{A}$ and
$\kappa_{L}$. Then the ratio between the maximum and minimum (non-zero)
singular value of $A_{L}$ is upper bounded by

$\kappa=\kappa_{L}\left(1+\frac{\norm{A}}{\sqrt{\lambda}\norm{L}}\right)$

We can also bound the spectral norm as

$\norm{A_{L}}=\Theta\left(\norm{A}+\sqrt{\lambda}\norm{L}\right)$

###### Proof.

To bound the spectral norm and condition number of $A_{L}$, consider the
eigenvalues of the following matrix:

$A_{L}^{T}A_{L}=\begin{pmatrix}A^{T}A+\lambda L^{T}L&0\\\ 0&0\end{pmatrix}$

This implies that the non-zero eigenvalues of $A_{L}^{T}A_{L}$ are the same as
those of $A^{T}A+\lambda L^{T}L$. Therefore, using triangle inequality, the
spectral norm of $A_{L}$ can be upper-bounded as follows:

$\displaystyle\norm{A_{L}}=\sqrt{\norm{A_{L}^{T}A_{L}}}=\sqrt{\norm{A^{T}A+\lambda
L^{T}L}}\leq\sqrt{\norm{A^{T}A}+\lambda\norm{L^{T}L}}=\sqrt{\norm{A}^{2}+\lambda\norm{L}^{2}}\leq\norm{A}+\sqrt{\lambda}\norm{L}$

Similarly $\norm{A_{L}}\geq\norm{A}$ and
$\norm{A_{L}}\geq\sqrt{\lambda}\norm{L}$, which effectively gives the tight
bound for $\norm{A_{L}}$.

As $L^{T}L$ is positive definite, we have that its minimum singular value is
$\sigma_{\mathrm{min}(L)}=\norm{L}/\kappa_{L}$. And we also know that $A^{T}A$
is positive semidefinite, so by Weyl’s inequality, the minimum singular value
of $A_{L}$ is lower bounded by

$\displaystyle\sigma_{\min}\left(A_{L}\right)\geq\sqrt{\sigma_{\min}\left(A\right)^{2}+\lambda\sigma_{\min}\left(L\right)^{2}}\geq\sqrt{\lambda\frac{\norm{L}^{2}}{\kappa_{L}^{2}}}=\sqrt{\lambda}\frac{\norm{L}}{\kappa_{L}}$

Thus,

$\dfrac{\sigma_{\max}\left(A_{L}\right)}{\sigma_{\min}\left(A_{L}\right)}\leq\kappa=\kappa_{L}\left(1+\frac{\norm{A}}{\sqrt{\lambda}\norm{L}}\right)$

∎

In the theorems and lemmas for regularized quantum linear regression and its
variants that we develop in this section, we consider that $L$ is a good
regularizer in order to provide a simple expression for $\kappa$. However,
this is without loss of generality. When $L$ is not a good regularizer, the
expressions for the respective complexities will remain unaltered, except that
$\kappa$ would now correspond to the condition number of the augmented matrix.

Now it might be possible that $\ket{b}$ does not belong to the row space of
$(A^{T}A+\lambda L^{T}L)^{-1}A^{T}$ which is equivalent to saying
$\ket{b}\ket{0}$ may not lie in $\mathrm{row}(A^{+}_{L})$. However, it is
reasonable to expect that the initial hypothesis of the underlying model being
close to linear is correct. That is, we expect $\ket{b}$ to have a good
overlap with
$\mathrm{row}\left(A^{+}_{L}\right)=\mathrm{col}\left(A_{L}\right)$. The
quantity that quantifies how far the model is from being linear is the so
called normalized residual sum of squares. For $\ell_{2}$-regularized ordinary
least squares, this is given by

$\mathcal{S}_{O}=\dfrac{\norm{(I-\Pi_{\mathrm{col}(A_{L})})\ket{b}\ket{0}}^{2}}{\norm{\ket{b}}^{2}}=1-\norm{\Pi_{\mathrm{col}(A_{L})}\ket{b}\ket{0}}^{2}.$
(33)

If the underlying data can indeed be fit by a linear function,
$\mathcal{S}_{O}$ will be low. Subsequently, we assume that
$\mathcal{S}_{O}=1-\norm{\Pi_{\mathrm{col}(A_{L})}\ket{b}\ket{0}}^{2}\leq\gamma<1/2$.
This in turn implies that
$\norm{\Pi_{\mathrm{col}(A_{L})}\ket{b}\ket{0}}^{2}=\Omega(1)$, implying that
the data can be reasonably fit by a linear model.iiiOur results also hold if
we assume that $\mathcal{S}_{O}\leq\gamma$ for some $\gamma\in(0,1)$. That is,
$\norm{\Pi_{\mathrm{col}(A_{L})}}\geq 1-\gamma$. In such a scenario our
complexity to prepare $A^{+}_{L}\ket{b,0}/\norm{A^{+}_{L}\ket{b,0}}$ is
rescaled by $1/\sqrt{1-\gamma}$.

Now we are in a position to present our quantum algorithm for the quantum
least squares problem with general $\ell_{2}$-regularization. We also present
an improved quantum algorithm for the closely related quantum ridge
regression, which is a special case of the former.

######  Theorem 32 (Quantum Ordinary Least Squares with General
$\ell_{2}$-Regularization).

Let $A,L\in\mathbb{R}^{N\times d}$ be the data and penalty matrices with
effective condition numbers $\kappa_{A}$ and $\kappa_{L}$ respectively, and
$\lambda\in\mathbb{R}^{+}$ be the regression parameter. Let $U_{A}$ be a
$(\alpha_{A},a_{A},\varepsilon_{A})$-block-encoding of $A$ implemented in time
$T_{A}$ and $U_{L}$ be a $(\alpha_{L},a_{L},\varepsilon_{L})$-block-encoding
of $L$ implemented in time $T_{L}$. Furthermore, suppose $U_{b}$ be a unitary
that prepares $\ket{b}$ in time $T_{b}$ and define

$\kappa=\kappa_{L}\left(1+\frac{\norm{A}}{\sqrt{\lambda}\norm{L}}\right)$

Then for any $\delta\in(0,1)$ such that

$\varepsilon_{A},\sqrt{\lambda}\varepsilon_{L}=o\left(\frac{\delta}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right)$
(34)

we can prepare a state that is $\delta$-close to

$\frac{\left(A^{T}A+\lambda
L^{T}L\right)^{-1}A^{T}\ket{b}}{\norm{\left(A^{T}A+\lambda
L^{T}L\right)^{-1}A^{T}\ket{b}}}$

with probability $\Theta(1)$, at a cost of

$\order{\kappa\log\kappa\left(\left(\frac{\alpha_{A}+\sqrt{\lambda}\alpha_{L}}{\norm{A}+\sqrt{\lambda}\norm{L}}\right)\log\left(\frac{\kappa}{\delta}\right)\left(T_{A}+T_{L}\right)+T_{b}\right)}$
(35)

using only $\order{\log\kappa}$ additional qubits.

###### Proof.

We invoke Lemma 22, to obtain a unitary $U$, which is a
$(\alpha_{A}+\sqrt{\lambda}\alpha_{L},\max(a_{A},a_{L})+2,\varepsilon_{A}+\sqrt{\lambda}\varepsilon_{L})$-block-
encoding of the matrix $A_{L}$, implemented at a cost of
$\order{T_{A}+T_{L}}$. Note that in Lemma 22, $A$ and $L$ are considered to be
$s$-qubit operators. For $N\times d$ matrices, such that $N,d\leq 2^{s}$, we
can pad them with zero entries. Padding $A$ and $L$ with zeros may result in
the augmented matrix $A_{L}$ having some zero rows between $A$ and $L$.
However, this is also not an issue as we are only interested in the top left
block of $A^{+}_{L}$ which remains unaffected.

Note that $U$ can be reinterpreted as a
$\left(\frac{\alpha_{A}+\sqrt{\lambda}\alpha_{L}}{\norm{A_{L}}},\max(a_{A},a_{L})+2,\frac{\varepsilon_{A}+\sqrt{\lambda}\varepsilon_{L}}{\norm{A_{L}}}\right)$-block-
encoding of the normalized matrix $A_{L}/\norm{A_{L}}$. Furthermore, we can
prepare the quantum state $\ket{b}\ket{0}$ in time $T_{b}$. Now by using
Theorem 28 with $U$ and an appropriately chosen $\delta$ specified above, we
obtain a quantum state that is $\delta$-close to

$\frac{(A^{T}A+\lambda L^{T}L)^{-1}A^{T}\ket{b}}{\norm{(A^{T}A+\lambda
L^{T}L)^{-1}A^{T}\ket{b}}}$

in the first register. ∎

In the above complexity, when $L$ is a good regularizer, $\kappa$ is
independent of $\kappa_{A}$. $\kappa$ can be made arbitrarily smaller than
$\kappa_{A}$ by an appropriate choice of $L$. Thus the regularized version has
significantly better time complexity than the unregularized case. One such
example of a good regularizer is in case of Quantum Ridge Regression, where we
use the identity matrix to regularize. The corollary below elucidates this.

######  Corollary 33 (Quantum Ridge Regression).

Let $A$ be a matrix of dimension $N\times d$ with effective condition number
$\kappa_{A}$ and $\lambda\in\mathbb{R}^{+}$ be the regression parameter. Let
$U_{A}$ be a $(\alpha,a,\varepsilon)$-block-encoding of $A$ implemented in
time $T_{A}$. Let $U_{b}$ be a unitary that prepares $\ket{b}$ in time
$T_{b}$. If $\kappa=1+\norm{A}/\sqrt{\lambda}$ then for any $\delta$ such that

$\varepsilon=o\left(\frac{\delta}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right)$

we can prepare a state $\delta$-close to

$\frac{\left(A^{T}A+\lambda
I\right)^{-1}A^{T}\ket{b}}{\norm{\left(A^{T}A+\lambda
I\right)^{-1}A^{T}\ket{b}}}$

at a cost of

$\order{\log\kappa\left(\frac{\alpha_{A}}{\sqrt{\lambda}}\log\left(\frac{\kappa}{\delta}\right)T_{A}+\kappa
T_{b}\right)}$ (36)

with probability $\Theta(1)$ using only $\order{\log\kappa}$ additional
qubits.

###### Proof.

The identity matrix $I$ is a trivial $(1,0,0)$-block-encoding of itself, and
$\kappa_{I}=1$. We invoke Theorem 32 with $L=I$ to obtain the solution. ∎

Being in the block-encoding framework allows us to express the complexity of
our quantum algorithm in specific input models such as the quantum data
structure input model and the sparse access model. We express these
complexities via the following corollaries.

######  Corollary 34 (Quantum Ordinary Least Squares with
$\ell_{2}$-Regularization in the Quantum Data Structure Input Model).

Let $A,L\in\mathbb{R}^{N\times d}$ with effective condition numbers
$\kappa_{A},\kappa_{L}$ respectively. Let $\lambda\in\mathbb{R}^{+}$ and
$b\in\mathbb{R}^{N}$. Let $\kappa$ be the effective condition number of the
augmented matrix $A_{L}$. Suppose that $A$, $L$ and $b$ are stored in a
quantum accessible data structure. Then for any $\delta>0$ there exists a
quantum algorithm to prepare a quantum state $\delta$-close to

$\frac{\left(A^{T}A+\lambda
L^{T}L\right)^{-1}A^{T}\ket{b}}{\norm{\left(A^{T}A+\lambda
L^{T}L\right)^{-1}A^{T}\ket{b}}}{}$

with probability $\Theta(1)$, at a cost of

$\order{\kappa\left(\frac{\mu_{A}+\sqrt{\lambda}\mu_{L}}{\norm{A}+\sqrt{\lambda}\norm{L}}\right)\mathrm{polylog}\left(Nd,\kappa,\frac{1}{\delta},\lambda\right)}.$
(37)

###### Proof.

Since $b$ is stored in the data structure, for some $\varepsilon_{b}>0$, we
can prepare the state $\ket{b^{\prime}}$ that is $\varepsilon_{b}$-close to
$\ket{b}=\sum_{i}b_{i}\ket{i}/\norm{b}$ using
$T_{b}=\order{\mathrm{polylog}(N/\varepsilon_{b})}$ queries to the data
structure (see Section 2.2.1.) Similarly, for some parameters
$\varepsilon_{A},\varepsilon_{L}>0$, we can construct a
$\left(\mu_{A},\lceil\log(d+N)\rceil,\varepsilon_{A}\right)$-block-encoding of
$A$ using $T_{A}=\order{\mathrm{polylog}(Nd/\varepsilon_{A})}$ queries to the
data structure and a
$\left(\mu_{L},\lceil\log(d+N)\rceil,\varepsilon_{B}\right)$-block-encoding of
$L$ using $T_{L}=\order{\mathrm{polylog}(Nd/\varepsilon_{B})}$ queries.

We invoke Theorem 32 with a precision $\delta/2$ by choosing $\varepsilon_{A}$
and $\varepsilon_{L}$ such that equation Equation 34 is satisfied. This gives
us a state that is $\delta/2$-close to

$\frac{\left(A^{T}A+\lambda
L^{T}L\right)^{-1}A^{T}\ket{b^{\prime}}}{\norm{\left(A^{T}A+\lambda
L^{T}L\right)^{-1}A^{T}\ket{b^{\prime}}}}$

To compute the final precision as $\delta$, we use Lemma 15 by choosing
$\varepsilon_{b}=\frac{\delta}{2\kappa}$. The complexity can be calculated by
plugging in the relevant values in Equation 35 ∎

In the previous corollary $\mu_{A}=\norm{A}_{F}$ and $\mu_{L}=\norm{L}_{F}$
when the matrix $A$ and $L$ are stored in the data structure. Similarly,
$\mu_{A}=\mu_{p}(A)$ and $\mu_{L}=\mu_{p}(L)$ when the matrices
$A^{(p)},A^{(1-p)}$ and $L^{(p)},L^{(1-p)}$ are stored in the data structure.

Now we discuss the complexity of quantum ordinary least squares with
$\ell_{2}$-regularization in the sparse access input model. We call a matrix
$M$ as $(s_{r},s_{c})$ row-column sparse if it has a row sparsity $s_{r}$ and
column sparsity $s_{c}$.

######  Corollary 35 (Quantum Ordinary least squares with
$\ell_{2}$-regularization in the sparse access model).

Let $A\in\mathbb{R}^{N\times d}$ be $(s^{A}_{r},s^{A}_{c})$ row-column sparse,
and similarly, let $L\in\mathbb{R}^{N\times d}$ be $(s^{L}_{r},s^{L}_{c})$
row-column sparse, with effective condition numbers $\kappa_{A}$ and
$\kappa_{L}$ respectively. Let $\lambda\in\mathbb{R}^{+}$ and $\delta>0$.
Suppose there exists a unitary that prepares $\ket{b}$ at a cost, $T_{b}$.
Then there is a quantum algorithm to prepare a quantum state that is
$\delta$-close to

$\frac{(A^{T}A+\lambda L^{T}L)^{-1}A^{T}\ket{b}}{\norm{(A^{T}A+\lambda
L^{T}L)^{-1}A^{T}\ket{b}}}$

with probability $\Theta(1)$, at a cost of

$\order{\kappa\left(\frac{\sqrt{s^{A}_{r}s^{A}_{c}}+\sqrt{\lambda
s^{L}_{r}s^{L}_{c}}}{\norm{A}+\sqrt{\lambda}\norm{L}}\right)\mathrm{polylog}\left(Nd,\kappa,\frac{1}{\delta},\lambda\right)+\kappa\log\kappa
T_{b}}.$ (38)

###### Proof.

The proof is similar to 34 but with $\alpha_{A}=\sqrt{s^{A}_{r}s^{A}_{c}}$ and
$\alpha_{L}=\sqrt{s^{L}_{r}s^{L}_{c}}$. ∎

###  4.2 Quantum Weighted And Generalized Least Squares

This technique of working with a augmented matrix will also hold for the other
variants of ordinary least squares. In this section, we begin by briefly
describing these variants before moving on to designing quantum algorithms for
the corresponding problems.

Weighted Least Squares:  For the WLS problem, each observation
$\\{a_{i},b_{i}\\}$ is assigned some weight $w_{i}\in\mathbb{R}^{+}$ and the
objective function to be minimized is of the form

$\mathcal{L}_{W}:=\sum_{j}w_{j}(x^{T}a_{j}-b_{j})^{2}.$ (39)

If $W\in\mathbb{R}^{N\times N}$ is the diagonal matrix with $w_{i}$ being the
$i^{\mathrm{th}}$ diagonal entry, then the optimal $x$ satisfies

$x=(A^{T}WA)^{-1}A^{T}Wb.$ (40)

The $\ell_{2}$-regularized version of WLS satisfies

$x=(A^{T}WA+\lambda L^{T}L)^{-1}A^{T}Wb$ (41)

Our quantum algorithm outputs a state that is close to

$\ket{x}=\frac{(A^{T}WA+\lambda
L^{T}L)^{-1}A^{T}W\ket{b}}{\norm{(A^{T}WA+\lambda L^{T}L)^{-1}A^{T}W\ket{b}}}$
(42)

given approximate block-encodings of $A$, $W$ and $L$. Much like Equation 32,
finding the optimal solution reduces to finding the pseudo inverse of an
augmented matrix $A_{L}$ given by

$A_{L}:=\begin{pmatrix}\sqrt{W}A&0\\\ \sqrt{\lambda}L&0\end{pmatrix}.$

The top left block of $A^{+}_{L}=(A^{T}WA+\lambda L^{T}L)^{-1}A^{T}\sqrt{W}$,
which is the required linear transformation to be applied to the vector
$y=\sqrt{W}b$. The ratio between the minimum and maximum singular values of
$A_{L}$, $\kappa$, can be obtained analogously to Lemma 31. For the
$\ell_{2}$-regularized WLS problem, normalized residual sum of squares is
given by

$\mathcal{S}_{W}=\dfrac{\norm{(I-\Pi_{\mathrm{col}(A_{L})})\ket{y}\ket{0}}^{2}}{\norm{\ket{y}}^{2}}=1-\norm{\Pi_{\mathrm{col}(A_{L})}\ket{y}\ket{0}}^{2}.$
(43)

Subsequently, we assume that
$\mathcal{S}_{W}=1-\norm{\Pi_{\mathrm{col}(A_{L})}\ket{y}\ket{0}}^{2}\leq\gamma<1/2$.
This in turn implies that
$\norm{\Pi_{\mathrm{col}(A_{L})}\ket{y}\ket{0}}^{2}=\Omega(1)$, implying that
the data can be reasonably fit by a linear model.

Generalized Least Squares. Similarly, we can extend this to GLS problem, where
there the input data may be correlated. These correlations are given by the
non-singular covariance matrix $\Omega\in\mathbb{R}^{N\times N}$. The WLS
problem is a special case of the GLS problem, corresponding to when $\Omega$
is a diagonal matrix. The objective function to be minimized is

$\mathcal{L}_{\Omega}:=\sum_{i,j}(\Omega^{-1})_{ij}(x^{T}a_{i}-b_{i})(x^{T}a_{j}-b_{j}).$
(44)

The optimal $x\in\mathbb{R}^{d}$ satisfies

$x=(A^{T}\Omega^{-1}A)^{-1}A^{T}\Omega^{-1}b$ (45)

Similarly, the $\ell_{2}$-regularized GLS solver outputs $x$ such that

$x=(A^{T}\Omega^{-1}A+\lambda L^{T}L)^{-1}A^{T}\Omega^{-1}b.$ (46)

So, given approximate block-encodings of $A$, $\Omega$ and $L$ a quantum GLS
solver outputs a quantum state close to

$\ket{x}=\frac{(A^{T}\Omega^{-1}A+\lambda
L^{T}L)^{-1}A^{T}\Omega^{-1}\ket{b}}{\norm{(A^{T}\Omega^{-1}A+\lambda
L^{T}L)^{-1}A^{T}\Omega^{-1}\ket{b}}}$ (47)

The augmented matrix $A_{L}$ is defined as

$A_{L}:=\begin{pmatrix}\Omega^{-1/2}A&0\\\ \sqrt{\lambda}L&0\end{pmatrix}.$

Then top left block of $A^{+}_{L}$ to the vector $y=\Omega^{-1/2}b$ yields the
optimal $x$. Thus the quantum GLS problem with $\ell_{2}$-regularization first
prepares $\Omega^{-1/2}\ket{b}\ket{0}$ and then uses the matrix inversion
algorithm by QSVT to implement $A^{+}_{L}\Omega^{-1/2}\ket{b}\ket{0}$.
Analogous to OLS and WLS, we assume that the normalized residual sum of
squares $\mathcal{S}_{\Omega}\leq\gamma<1/2$.

####  4.2.1 Quantum Weighted Least Squares

In this section, we derive the complexity of the $\ell_{2}$-regularized WLS
problem. We assume that we have a diagonal weight matrix
$W\in\mathbb{R}^{N\times N}$ such that its smallest and largest diagonal
entries are ${w_{\text{min}}}$ and ${w_{\text{max}}}$, respectively. This
implies that $\norm{W}={w_{\text{max}}}$ and
$\kappa_{W}={w_{\text{max}}}/{w_{\text{min}}}$. We take advantage of the fact
that the matrix $W$ is diagonal and then apply controlled rotations to
directly implement a block encoding of $\sqrt{W}A$. Additionally, given a
state preparation procedure for $\ket{b}$, we can easily prepare a state
proportional to $\sqrt{W}\ket{b}$. We then use Theorem 32 to solve QWLS.

We first formalize this idea in Theorem 36, assuming direct access to (i) a
block encoding of $B=\sqrt{W}A$, and (ii) a procedure for preparing the state
$\ket{b_{w}}=\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$. Subsequently,
for the specific input models, we show that we can indeed efficiently obtain a
block-encoding of $B$ and prepare the state $\ket{b_{w}}$.

######  Theorem 36 (Quantum Weighted Least Squares with General
$\ell_{2}$-Regularization).

Let $A,L\in\mathbb{R}^{N\times d}$, be the data and penalty matrix, with
effective condition numbers $\kappa_{A}$ and $\kappa_{L}$, respectively. Let
$\lambda\in\mathbb{R}^{+}$ be the regularizing parameter. Let
$W\in\mathbb{R}^{N\times N}$ be a diagonal weight matrix with the largest and
smallest diagonal entries being ${w_{\text{max}}},{w_{\text{min}}}$,
respectively. Let $U_{B}$ be a $(\alpha_{B},a_{B},\varepsilon_{B})$ block
encoding of $B:=\sqrt{W}A$ implemented in time $T_{B}$ and let $U_{L}$ be a
$(\alpha_{L},a_{L},\varepsilon_{L})$ block encoding of $L$ implemented in time
$T_{L}$, such that
$\varepsilon_{B}=o\left(\frac{\delta}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right)$
and
$\varepsilon_{L}=o\left(\frac{\delta}{\sqrt{\lambda}\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right)$.
Let $U_{b_{w}}$ be a unitary that prepares
$\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$ in time $T_{b_{w}}$. Define

$\kappa:=\kappa_{L}\left(1+\frac{\sqrt{{w_{\text{max}}}}\norm{A}}{\sqrt{\lambda}\norm{L}}\right)$

Then for any $\delta>0$ we can prepare a quantum state that is $\delta$-close
to

$\frac{(A^{T}WA+\lambda L^{T}L)^{-1}A^{T}W\ket{b}}{\norm{(A^{T}WA+\lambda
L^{T}L)^{-1}A^{T}W\ket{b}}}$

with probability $\Theta(1)$, at a cost of

$\order{\kappa\log\kappa\left(\frac{\alpha_{B}+\sqrt{\lambda}\alpha_{L}}{\sqrt{{w_{\text{max}}}}\norm{A}+\sqrt{\lambda}\norm{L}}\log\left(\frac{\kappa}{\delta}\right)(T_{B}+T_{L})+T_{b_{w}}\right)},$
(48)

using only $\order{\log\kappa}$ additional qubits.

###### Proof.

We then invoke Theorem 32 with $B$ and $L$ as the data and regularization
matrices, respectively. This requires that $\varepsilon_{B},\varepsilon_{L}$
such that

$\varepsilon_{B}+\sqrt{\lambda}\varepsilon_{L}=o\left(\frac{\delta}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta}\right)}\right).$

Thus, we get the upper bounds on the precision
$\varepsilon_{B},\leavevmode\nobreak\ \varepsilon_{L}$ required. This gives us
a quantum state $\delta$-close to

$\frac{(A^{T}WA+\lambda L^{T}L)^{-1}A^{T}W\ket{b}}{\norm{(A^{T}WA+\lambda
L^{T}L)^{-1}A^{T}W\ket{b}}}.$

∎

Next, we construct the block encodings for $\sqrt{W}A$ and the state
$\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$ efficiently in the quantum
data structure input model. This construction would also apply to the sparse
access input model with slight modifications.

######  Lemma 37 (Efficiently preparing $\sqrt{W}A$ in the Quantum Data
Structure Model).

Let $W\in\mathbb{R}^{N\times N}$ such that $W=\mathrm{diag}(w_{1},w_{2}\ldots
w_{N})$ and $w_{\max}:=\max_{i}w_{i}$, and $A\in\mathbb{R}^{N\times d}$ be
stored in a quantum-accessible data structure. Then for any $\delta>0$ there
exists a

$(\sqrt{w_{\max}}\norm{A}_{F},\lceil\log\left(N+d\right)\rceil,\delta)$

block-encoding of $\sqrt{W}A$ that can be implemented at the cost
$\order{\mathrm{polylog}(Nd/\delta)}$.

###### Proof.

$\forall j\in[N]$, define

$\ket{\psi_{j}}:=\sqrt{\frac{w_{j}}{w_{\max}}}\ket{j}\frac{1}{\norm{A_{j,\cdot}}}\sum_{k\in[d]}A_{j,k}\ket{k}.$

Similarly, $\forall k\in[d]$, define

$\ket{\phi_{k}}:=\frac{1}{\norm{A}_{F}}\left(\sum_{j\in[N]}\norm{A_{j,\cdot}}\ket{j}\right)\ket{k}.$

Observe that $\forall j\in[N],k\in[d]$,

$\braket{\psi_{j}}{\phi_{k}}=\sqrt{\frac{w_{j}}{w_{\max}}}\frac{A_{j,k}}{\norm{A}_{F}}=\frac{\bra{j}\sqrt{W}A\ket{k}}{\sqrt{w_{\max}}\norm{A}_{F}}.$

Given quantum data structure accesses to $W$ and $A$, one can construct
quantum circuits $W_{R}$ and $W_{L}$ similar to $U_{L}$ and $U_{R}$ from Lemma
3 that prepare $\ket{\phi_{k}}$ and $\ket{\psi_{j}}$ above. $\ket{\phi_{k}}$
can be prepared just as in Lemma 3, while $\ket{\psi_{j}}$ can be prepared
using controlled rotations on the state $\ket{\frac{w_{j}}{w_{\max}}}$ (which
can be constructed from the QRAM access to $W$) after adding an ancilla qubit
and the QRAM access to $A$. Thus, $W_{R}^{\dagger}W_{L}$ is the required block
encoding, which according to Theorem 2 can be implemented using
$\mathrm{polylog}(Nd/\delta)$ queries. ∎

######  Lemma 38 (Efficiently preparing $\sqrt{W}\ket{b}$ in the Quantum Data
Structure Model).

Let $b\in\mathbb{R}^{N}$ and $W\in\mathbb{R}^{N\times N}$. Suppose that $b$
and $W$ are stored in a quantum-accessible data structure such that we have a
state preparation procedure that acts as

$\displaystyle U_{W}:\ket{j}\ket{0}\mapsto\ket{j}\ket{w_{j}},$ $\displaystyle
U_{b}:\ket{0}\mapsto\sum_{j}\frac{b_{j}}{\norm{b}}\ket{j}.$

Then for any $\delta>0$ we can prepare the quantum state that is
$\delta$-close to $\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$ with
constant success probability and at a cost of
$\order{\sqrt{\frac{w_{\max}}{w_{\min}}}\mathrm{polylog}\left(\frac{N}{\delta}\right)}$.

###### Proof.

Use $U_{b}$ to prepare the state

$\ket{b}=\frac{1}{\norm{b}}\sum_{j}b_{j}\ket{j}$

in time polylog$(N)$. Then, apply the following transformation

$\displaystyle\ket{j}\ket{0}\ket{0}$
$\displaystyle\mapsto\ket{j}\ket{w_{j}}\ket{0}$
$\displaystyle\mapsto\ket{j}\ket{w_{j}}\left(\sqrt{\frac{w_{j}}{w_{\max}}}\ket{0}+\sqrt{1-\frac{w_{j}}{w_{\max}}}\ket{1}\right)$
$\displaystyle\mapsto\ket{j}\ket{0}\left(\sqrt{\frac{w_{j}}{w_{\max}}}\ket{0}+\sqrt{1-\frac{w_{j}}{w_{\max}}}\ket{1}\right)$

which can again be applied using some controlled rotations, a square root
circuit and $U_{W}$. This gives us the state (ignoring some blank registers)

$\sum_{j}\left(\sqrt{\frac{w_{j}}{w_{\max}}}\ket{0}+\sqrt{1-\frac{w_{j}}{w_{\max}}}\ket{1}\right)\frac{b_{j}}{\norm{b}}\ket{j}.$
(49)

The probability for the ancilla to be in $\ket{0}$ state is

$\Omega\left(\frac{w_{\min}}{w_{\max}}\right).$

Thus performing $\order{\sqrt{\frac{w_{\max}}{w_{\min}}}}$ rounds of amplitude
amplification on $\ket{0}$ gives us a constant probability of observing
$\ket{0}$, and therefore obtaining the desired state
$\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$. ∎

Using the above two theorems, and the quantum OLS solver (Theorem 32), we can
construct an algorithm for regularized quantum WLS.

######  Theorem 39 (Quantum Weighted Least Squares with General
$\ell_{2}$-Regularization in the Quantum Data Structure Model).

Let $A,L\in\mathbb{R}^{N\times d}$ with effective condition numbers
$\kappa_{A},\kappa_{L}$ respectively be stored in an efficient quantum
accessible data structure. Let $W\in\mathbb{R}^{N\times N}$ be a diagonal
matrix with largest and smallest singular values
${w_{\text{max}}},{w_{\text{min}}}$ respectively, which is also stored in an
efficient quantum accessible data structure. Furthermore, suppose the entries
of the vector $b\in\mathbb{R}^{N}$ are also stored in a quantum-accessible
data structure and define,

$\kappa:=\kappa_{L}\left(1+\frac{\sqrt{{w_{\text{max}}}}\norm{A}}{\sqrt{\lambda}\norm{L}}\right)$

Then for any $\delta>0$ we can prepare a quantum state that is $\delta$-close
to

$\frac{(A^{T}WA+\lambda L^{T}L)^{-1}A^{T}W\ket{b}}{\norm{(A^{T}WA+\lambda
L^{T}L)^{-1}A^{T}W\ket{b}}}$

with probability $\Theta(1)$, at a cost of

$\order{\kappa\left(\frac{\sqrt{w_{\text{max}}}\norm{A}_{F}+\sqrt{\lambda}\norm{L}_{F}}{\sqrt{w_{\text{max}}}\norm{A}+\sqrt{\lambda}\norm{L}}+\sqrt{\frac{{w_{\text{max}}}}{{w_{\text{min}}}}}\right)\mathrm{polylog}\left(Nd,\kappa,\frac{1}{\delta}\right)}$
(50)

###### Proof.

Choose some precision parameter $\varepsilon>0$ for accessing the data
structure. Given access to $W$ and $A$, we can use Lemma 37 to prepare a
$(\sqrt{{w_{\text{max}}}}\norm{A}_{F},\lceil\log\left(N+d\right)\rceil,\varepsilon)$-block-
encoding of $\sqrt{W}A$, using
$T_{A}:=\order{\mathrm{polylog}\left(Nd/\varepsilon\right)}$ queries to the
data structure. Similarly, Lemma 3 allows us to build a
$(\norm{L}_{F},\lceil\log\left(N+d\right)\rceil,\varepsilon)$-block-encoding
of $L$ using $T_{L}:=\order{\mathrm{polylog}(Nd/\varepsilon)}$ queries to the
data structure.

Next, using Lemma 38, for any $\varepsilon_{b}>0$, we can prepare a state
$\varepsilon_{b}$-close to
$\ket{b^{\prime}}:=\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$. This
procedure requires
$T_{b}:=\order{\sqrt{\frac{w_{\max}}{w_{\min}}}\mathrm{polylog}\left(N/\varepsilon_{b}\right)}$
queries to the data structure. Now we can invoke the OLS solver in Theorem 32
with a precision of $\delta_{b}$, by considering $\sqrt{W}A$ as the data
matrix and $\frac{\sqrt{W}\ket{b}}{\norm{\sqrt{W}\ket{b}}}$ as the input
state. In order for the input block-encoding precision to satisfy the bound in
Equation 34, we choose $\varepsilon$ such that

$\varepsilon=o\left(\frac{\delta_{b}}{\kappa^{3}\log^{2}\left(\frac{\kappa}{\delta_{b}}\right)}\right).$

Finally, for the output state to be $\delta$-close to the required state, we
choose $\delta_{b}=\delta/2$ and $\varepsilon_{b}=\delta/2\kappa$ to use the
robustness result from Lemma 15. This gives us