# Quantum circuit design methodology for multiple linear regression

Sanchayan Dutta, Adrien Suau, Sagnik Dutta, Suvadeep Roy, Bikash K. Behera, Prasanta K. Panigrahi

► To cite this version:

Sanchayan Dutta, Adrien Suau, Sagnik Dutta, Suvadeep Roy, Bikash K. Behera, et al., Quantum circuit design methodology for multiple linear regression. IET Quantum Communication, 2020, 1 (2), pp. 55-61. 10.1049/iet-qtc.2020.0013. lirmm-03035852

HAL Id: lirmm-03035852

https://hal-lirmm.ccsd.cnrs.fr/lirmm-03035852

Submitted on 2 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantum circuit design methodology for multiple linear regression

eISSN 2632-8925
Received on 21st July 2020
Revised 26th September 2020
Accepted on 7th October 2020
doi: 10.1049/iet-qtc.2020.0013
www.ietdl.org

Sanchayan Dutta<fnref n="1" />, Adrien Suau<fnref n="2" />, Sagnik Dutta<fnref n="3" />, Suvadeep Roy<fnref n="3" />, Bikash K. Behera<fnref n="3" />,<fnref n="4" />, Prasanta K. Panigrahi<fnref n="3" />

**Abstract:** Multiple linear regression assumes an imperative role in supervised machine learning. In 2009, Harrow et al. [Phys. Rev. Lett. 103, 150502 (2009)] showed that their Harrow Hassidim Lloyd (HHL) algorithm can be used to sample the solution of a linear system $Ax = b$ exponentially faster than any existing classical algorithm. The entire field of quantum machine learning gained considerable traction after the discovery of this celebrated algorithm. However, effective practical applications and experimental implementations of HHL are still sparse in the literature. Here, the authors demonstrate a potential practical utility of HHL, in the context of regression analysis, using the remarkable fact that there exists a natural reduction of any multiple linear regression problem to an equivalent linear systems problem. They put forward a 7-qubit quantum circuit design, motivated from an earlier work by Cao et al. [Mol. Phys. 110, 1675 (2012)], to solve a three-variable regression problem, using only elementary quantum gates. They also implement the group leaders optimisation algorithm (GLOA) [Mol. Phys. 109 (5), 761 (2011)] and elaborate on the advantages of using such stochastic algorithms in creating low-cost circuit approximations for the Hamiltonian simulation. Further, they discuss their Qiskit simulation and explore certain generalisations to the circuit design.

## 1 Introduction

Quantum algorithms running on quantum computers aim at quickly and efficiently solving several important computational problems faster than classical algorithms running on classical computers [1–11]. One key way in which quantum algorithms differ from classical algorithms is that they utilise quantum mechanical phenomena such as superposition and entanglement, that allows us to work in exponentially large Hilbert spaces with only polynomial overheads. This in turn, in some cases, allows for exponential speed-ups in terms of algorithmic complexity [1].

In today's world, machine learning is primarily concerned with the development of low-error models in order to make accurate predictions possible by learning and inferring from training data [12, 13]. It borrows heavily from the field of statistics in which linear regression is one of the flagship tools. The theory of multiple linear regression or more generally multivariate linear regression was largely developed in the field of statistics in the pre-computer era. It is one of the most well understood, versatile and straightforward techniques in any statistician's toolbox. It is also an important and practical supervised learning algorithm. Supervised learning is where one has some labelled input data samples $\{x_i, y_i\}_{i=1}^N$ (where $x_i$'s are the feature vectors and $y_i$'s are the corresponding labels) and then based on some criteria (which might depend on the context) chooses a mapping from the input set $X$ to the output set $Y$. That mapping can help to predict the probable output corresponding to an input lying outside of the training data sets. Multiple linear regression is similar in the sense that given some training samples one identifies a closely fitting hyperplane depending on the specific choice of a loss function (the most common one being a quadratic loss function based on the ‘least squares’ method). Interestingly, any multiple regression problem can be converted into an equivalent system of linear equations problem or more specifically, a quantum linear systems problem (QLSP) [14]. The process has been outlined using an example in Section 3.

Suppose that we are given a system of $N$ linear equations with $N$ unknowns, which can be expressed as $Ax = b$. Now, what we are interested in, is: given a matrix $A \in \mathbb{C}^{N \times N}$ with a vector $b \in \mathbb{C}^N$, find the solution $x \in \mathbb{C}^N$ satisfying $Ax = b$ (which is $A^{-1}b$ if $A$ is invertible), or else return a flag if no solution exists. This is known as the linear systems problem (LSP). However, we will consider only a special case of this general problem, in the form of the QLSP [14, 15].

The quantum version of the LSP problem, the QLSP, can be expressed as follows.

Let $A$ be a $N \times N$ Hermitian matrix with a spectral norm bounded by unity and a known condition number $\kappa$. The quantum state on $[\log N]$ qubits $|b\rangle$ can be given by

$$|b\rangle := \frac{\sum_i b_i |i\rangle}{\| \sum_i b_i |i\rangle \|} \qquad (1)$$

and $|x\rangle$ by

$$|x\rangle := \frac{\sum_i x_i |i\rangle}{\| \sum_i x_i |i\rangle \|} \qquad (2)$$

where $b_i, x_i$ are, respectively, the $i$th component of vectors $b$ and $x$. Given the matrix $A$ (whose elements are accessed by an oracle) and the state $|b\rangle$, an output state $|\tilde{x}\rangle$ is such that $\| |\tilde{x}\rangle - |x\rangle \|_2 \le \epsilon$, with some probability $\Omega(1)$ (practically at least 1/2) along with a binary flag indicating ‘success’ or ‘failure’ [14].

The restrictions on Hermiticity and spectral norm can be relaxed by noting that, even for a non-Hermitian matrix $A$, the corresponding

$$\begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix}$$

<footnote n="1">¹Department of Electronics and Telecommunication Engineering, Jadavpur University, Kolkata 700032, West Bengal, India</footnote>

<footnote n="2">²CERFACS, 42 Avenue Gaspard Coriolis, 31100 Toulouse, France</footnote>

<footnote n="3">³Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur 741246, West Bengal, India</footnote>

<footnote n="4">⁴Bikash's Quantum (OPC) Private Limited, Balindi, Mohanpur 741246, Nadia, West Bengal, India</footnote>

<footnote>✉ E-mail: sanchayan98@gmail.com</footnote>
matrix is Hermitian. This implies that we can instead solve the linear system given by

$$ \begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix} y = \begin{bmatrix} b \\ 0 \end{bmatrix}, $$

which has the unique solution

$$ y = \begin{bmatrix} 0 \\ x \end{bmatrix} $$

when **A** is invertible [4]. Also, any non-singular matrix can be scaled appropriately to adhere to the given conditions on the eigenspectrum. Note that the case when **A** is non-invertible has already been excluded by the fact that a known finite condition number exists for the matrix **A**.

In 2009, Harrow, Hassidim and Lloyd [4] put forward a quantum algorithm (popularly known as the 'HHL algorithm') to obtain information about the solution **x** of certain classes of linear systems **Ax** = **b**. As we know, algorithms for finding the solutions to linear systems of equations play an important role in engineering, physics, chemistry, computer science and economics apart from other areas. The HHL algorithm is highly celebrated in the world of quantum algorithms and quantum machine learning, but there is still a distinct lack of literature demonstrating practical uses of the algorithm and their experimental implementations. Experimentally implementing the HHL algorithm for solving an arbitrary system of linear equations to a satisfactory degree of accuracy remains an infeasible task even today, due to several physical and theoretical restrictions imposed by the algorithm and the currently accessible hardware. Our fundamental aim in this paper is to put forth a neat technique to demonstrate how HHL can be used for one of the most useful tools in statistical modelling – multiple regression analysis. We show that there is a canonical reduction of any multiple linear regression problem into a QLSP that is thereafter amenable to a solution using HHL. It turns out that, after slight modifications and classical pre-processing, the circuit used for HHL can be used for multiple linear regression as well. Here we should point out, in Scott Aaronson's words [16], that the HHL algorithm is mostly a template for other quantum algorithms and it has some caveats that must be kept in mind during the circuit design and physical implementation step, to preserve the exponential speedup offered by the algorithm. Nonetheless, these issues are not always detrimental and can often be overcome by appropriately preparing and justifying the circuit design methodology for specific purposes, as we will do here. There is another brief discussion related to these issues in the conclusions section of the paper.

In Section 4, we present an application (in the context of multiple regression) of a modified version of the earlier circuit design by Cao et al. [17] which was meant for implementing the HHL algorithm for a 4 × 4 linear system on real quantum computers. This circuit requires only 7 qubits and it should be simple enough to experimentally verify it if one gets access to a quantum processor having logic gates with sufficiently low error rates. Previously, Pan et al. demonstrated HHL on a 4-qubit NMR quantum computer [18], so we believe that it will be easily possible to experimentally implement the circuit we discuss, given the rapid rise in the number of qubits in quantum computer chips, in the past few years.

We also note that although HHL solves the QLSP for all matrices **A** or

$$ \begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix}, $$

it can be efficiently implemented only when they are sparse and well-conditioned (the sparsity condition may be slightly relaxed) [14]. In this context, 'efficient' means 'at most polylogarithmic in system size'. A $N \times N$ matrix is called *s-sparse* if it has at most *s* non-zero entries in any row or column. We call it simply *sparse* if it has at most poly(log *N*) entries per row [15]. We generally call a

matrix well-conditioned when its singular values lie between the reciprocal of its condition number (1/κ) and 1 [4]. Condition number κ of a matrix is the ratio of largest to smallest singular value and is undefined when the smallest singular value of **A** is 0. For Hermitian matrices, the eigenvalue magnitudes are equal to the magnitudes of the respective singular values.

At this point, it is important to reiterate that unlike the output $A^{-1}b$ of a classical linear system solver, the output copy of $|\bar{x}\rangle$ does not provide access to the coordinates of $A^{-1}b$. Nevertheless, it allows for sampling from the solution vectors like $\langle\bar{x}|M|\bar{x}\rangle$, where $M$ is a quantum-mechanical operator. This is one main difference between solutions of the LSP and solutions of the QLSP. We should also keep in mind that reading out the elements of $|\bar{x}\rangle$ in itself takes $\mathcal{O}(N)$ time. Thus, a solution to QLSP might be useful only in applications where just samples from the vector $|\bar{x}\rangle$ are needed [4, 14, 15].

The best existing classical matrix inversion algorithm involves the Gaussian elimination technique which takes $\mathcal{O}(N^3)$ time. For s-sparse and positive semi-definite **A**, the conjugate gradient algorithm [19] can be used to find the solution vector **x** in $\mathcal{O}(Nsklog(1/\epsilon))$ time by minimising the quadratic error function $|Ax - b|^2$, where *s* is the matrix sparsity, *κ* is the condition number and *ε* is the desired precision parameter. On the other hand, the HHL algorithm scales as $\mathcal{O}(\log(N)\epsilon^2\kappa^2/\epsilon)$, and is exponentially faster in *N* but polynomially slower in *s* and *κ*. In 2010, Andris Ambainis further improved the runtime of HHL to $\mathcal{O}(\kappa\log^3\kappa\log N/\epsilon)$ [20]. The exponentially worse slowdown in *ε* was also eliminated by Childs et al. in 2017 [14] and it got improved to $\mathcal{O}(skpolylog(sk/\epsilon))$ [15]. Since HHL has logarithmic scaling only for sparse or low-rank matrices, in 2018, Wossnig et al. [21] extended the HHL algorithm with quantum singular value estimation and provided a quantum linear system algorithm for dense matrices which achieves a polynomial improvement in time complexity, that is, $\mathcal{O}(\sqrt{N}\text{polylog}(N)\epsilon^2/\epsilon)$ (HHL retains its logarithmic scaling only for sparse or low-rank matrices). Furthermore, an exponential improvement is achievable with this algorithm if the rank of **A** is polylogarithmic in the matrix dimension.

Last but not least, we note that it is assumed that the state $|b\rangle$ can be efficiently constructed, i.e. prepared in 'polylogarithmic time'. In reality, however, efficient preparation of arbitrary quantum states is hard and is subject to several constraints.

## 2 Harrow Hassidim Lloyd algorithm

The HHL algorithm consists of three major steps which we will briefly discuss one by one. Initially, we begin with a Hermitian matrix **A** and an input state $|b\rangle$ corresponding to our specific system of linear equations. The assumption that **A** is Hermitian may be dropped without loss of generality since we can instead solve the linear system of equations given by

$$ \begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix} y = \begin{bmatrix} b \\ 0 \end{bmatrix} $$

which has the unique solution

$$ y = \begin{bmatrix} 0 \\ x \end{bmatrix} $$

when **A** is invertible. This transformation does not alter the condition number (ratio of the magnitudes of the largest and smallest eigenvalues) of **A** [4, 14]. However, in the case our original matrix **A** is not Hermitian, the transformed system with the new matrix

$$ \begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix} $$

needs oracle access to the non-zero entries of the rows and columns of **A** [14]. Since **A** is assumed to be Hermitian, it follows

The HHL algorithm consists of three major steps which we will briefly discuss one by one. Initially, we begin with a Hermitian matrix **A** and an input state $|b\rangle$ corresponding to our specific system of linear equations. The assumption that **A** is Hermitian may be dropped without loss of generality since we can instead solve the linear system of equations given by

$$ \begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix} y = \begin{bmatrix} b \\ 0 \end{bmatrix} $$

which has the unique solution

$$ y = \begin{bmatrix} 0 \\ x \end{bmatrix} $$

when **A** is invertible. This transformation does not alter the condition number (ratio of the magnitudes of the largest and smallest eigenvalues) of **A** [4, 14]. However, in the case our original matrix **A** is not Hermitian, the transformed system with the new matrix

$$ \begin{bmatrix} 0 & A^{\dagger} \\ A & 0 \end{bmatrix} $$
<figure><img src="image_2.png" /><figcaption>**Fig. 1** HHL algorithm schematic [The HHL algorithm schematic (Fig. 1) was generated using the TikZ code provided by Dr. Niel de Beaudrap (Department of Computer Science, Oxford University). It was inspired by figure 5 of the Dervovic et al. paper [15].]</figcaption></figure>

that $e^{iAt}$ is unitary. Here $iAt$ and $-iAt$ commute and hence $e^{iAt}e^{-iAt} = e^{iAt - iAt} = e^0 = I$. Moreover, $e^{iAt}$ shares all its eigenvectors with $A$, while its eigenvalues are $e^{i\lambda_j t}$ if the eigenvalues of $A$ are taken to be $\lambda_j$. Suppose that $|u_j\rangle$ are the eigenvectors of $A$ and $\lambda_j$ are the corresponding eigenvalues. We recall that we assumed all the eigenvalues to be of magnitude $<1$ (spectral norm is bounded by unity). As the eigenvalues $\tilde{\lambda}_j$ are of the form $0.a_1a_1a_3...$ in binary [1, p. 222], we will use $|\lambda_j\rangle$ to refer to $|a_1a_2a_3...\rangle$. We know from the spectral theorem that every Hermitian matrix has an orthonormal basis of eigenvectors. So, in this context, $A$ can be re-written as $\sum_j \lambda_j |u_j\rangle\langle u_j|$ (via eigendecomposition of $A$) and $|b\rangle$ as $\sum_j \beta_j |u_j\rangle$.

## 2.1 Phase estimation

The quantum phase estimation algorithm performs the mapping $((|0\rangle^{\otimes n})^C|u\rangle^I|0\rangle^S) \mapsto |\tilde{\phi}\rangle^C|u^I|0\rangle^S$ where $|u\rangle$ is an eigenvector of a unitary operator $U$ with an unknown eigenvalue $e^{i2\pi\phi}$ [1]. $\tilde{\phi}$ is a t-bit approximation of $\phi$, where $t$ is the number of qubits in the clock register. The superscripts on the kets indicate the names of the registers which store the corresponding states. In the HHL algorithm, the input register begins with a superposition of eigenvectors instead, i.e. $|b\rangle = \sum_j \beta_j |u_j\rangle$ instead of a specific eigenvector $|u\rangle$, and for us the unitary operator is $e^{iAt}$. So, the phase estimation circuit performs the mapping

$$((|0\rangle^{\otimes n})^C|b\rangle^I \mapsto \left( \sum_{j=1}^{N} \beta_j |u_j\rangle^I \left| \frac{\tilde{\lambda}_j f_0}{2\pi} \right) \right)^C)$$

where $\tilde{\lambda}_j$'s are the binary representations of the eigenvalues of $A$ to a tolerated precision. To be more explicit, here $\tilde{\lambda}_j$ is represented as $b_1b_2b_3...b_t$ (t being number of qubits in the clock register) if the actual binary equivalent of $\tilde{\lambda}_j$ is of the form $\lambda = 0.b_1b_2b_3...$. To avoid the factor of $2\pi$ in the denominator, the 'evolution time' $t_0$ is generally chosen to be $2\pi$. However, $t_0$ may also be used to 'normalise' $A$ (by re-scaling $t_0$) in case the spectral norm of $A$ exceeds 1 [Ideally, we should know both the upper bound and the lower bound of the eigenvalues, for effective re-scaling. Furthermore, to get accurate estimates, we should attempt to spread the possible values of $\lambda t$ over the whole $2\pi$ range.]. Additionally, an important factor in the performance of the algorithm is the condition number $\kappa$. As $\kappa$ grows, $A$ tends more and more towards a non-invertible matrix, and the solutions become less and less stable [4]. Matrices with large condition numbers are said to be 'ill-conditioned'. The HHL algorithm generally assumes that the singular values of $A$ lie between $1/\kappa$ and 1, which ensures that the matrices we have to deal with are 'well-conditioned'. Nonetheless, there are methods to tackle ill-conditioned matrices and those have been thoroughly discussed in the paper by Lloyd and co-workers [4]. It is worth mentioning that in this step the 'clock register'-controlled Hamiltonian simulation gate $U$ can be expressed as

$\sum_{k=0}^{T-1} |\tau\rangle\langle\tau|^C \otimes e^{iAt} t_0/T$, where $T=2^t$ and evolution time $t_0 = \mathcal{O}(\kappa/\epsilon)$. Interestingly choosing $t_0 = \mathcal{O}(\kappa/\epsilon)$ can at worse cause an error of magnitude $\epsilon$ in the final state [4].

## 2.2 R($\tilde{\lambda}^{-1}$) rotation

A 'clock register' controlled $\sigma_y$-rotation of the 'ancilla' qubit produces a normalised state of the form

$$\sum_{j=1}^{N} \beta_j |u_j\rangle^I |\tilde{\lambda}_j\rangle^C \left( \sqrt{1 - \frac{C^2}{\tilde{\lambda}_j^2}} |0\rangle + \frac{C}{\tilde{\lambda}_j} |1\rangle \right)^S$$

These rotations, conditioned on respective $\tilde{\lambda}_j$, can be achieved by the application of the

$$\exp(-i\theta\sigma_y) = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}$$

operators where $\theta = \cos^{-1}(C/\tilde{\lambda}_j)$. $C$ is a scaling factor to prevent the controlled rotation from being unphysical [15]. That is, practically $C < \lambda_{min}$ is a safe choice, which may be more formally stated as $C = \mathcal{O}(1/\kappa)$ [4] (Fig. 1).

## 2.3 Uncomputation

In the final step, the inverse quantum phase estimation algorithm sets back the clock register to $((|0\rangle^{\otimes n})^C)$ and leaves the remaining state as

$$\sum_{j=1}^{N} \beta_j |u_j\rangle^I \left( \sqrt{1 - \frac{C^2}{\tilde{\lambda}_j^2}} |0\rangle + \frac{C}{\tilde{\lambda}_j} |1\rangle \right)^S$$

Postselecting on the ancilla $|1\rangle^S$ gives the final state $C\sum_{j=1}^{N} (\beta_j/\tilde{\lambda}_j) |u_j\rangle^I$ [15]. The inverse of the Hermitian matrix $A$ can be written as $\sum_j (1/\tilde{\lambda}_j) |u_j\rangle\langle u_j|$, and hence $A^{-1}|b\rangle$ matches $\sum_j (\beta_j/\tilde{\lambda}_j) |u_j\rangle^I$. This outcome state, in the standard basis, is component-wise proportional to the exact solution $x$ of the system $Ax = b$ [17].

# 3 Linear regression utilising HHL

Linear regression models a linear relationship between a scalar 'response' variable and one or more 'feature' variables. Given a $n$-unit data set $\{y_i, x_{i1}, ..., x_{ip}\}_{i=1}^n$, a linear regression model assumes that the relationship between the dependent variable $y$ and a set of $p$ attributes, i.e. $x = \{x_1, ..., x_p\}$ is linear [12]. Essentially, the model takes the form

$$y_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip} + \epsilon_i = x_i^T \beta + \epsilon_i$$

Linear regression models a linear relationship between a scalar 'response' variable and one or more 'feature' variables. Given a $n$-unit data set $\{y_i, x_{i1}, ..., x_{ip}\}_{i=1}^n$, a linear regression model assumes that the relationship between the dependent variable $y$ and a set of $p$ attributes, i.e. $x = \{x_1, ..., x_p\}$ is linear [12]. Essentially, the model takes the form
where $\epsilon_i$ is the noise or error term. Here $i$ ranges from 1 to $n$. $x_i^T$ denotes the transpose of the column matrix $x_i$. $x_i^T \beta$ is the inner product between vectors $x_i$ and $\beta$. These $n$ equations may be more compactly represented in the matrix notation, as $y = X\beta + \epsilon$. Now, we will consider a simple example with three feature variables and a bias $\beta_0$. Say our data sets are
$$ \begin{aligned} & \{-(1/8) + (1/8\sqrt{2}), -\sqrt{2}, (1/\sqrt{2}), -(1/2)\}, \\ & \{(3/8) - (3/8\sqrt{2}), -\sqrt{2}, -(1/\sqrt{2}), (1/2)\}, \\ & \{-(1/8) - (1/8\sqrt{2}), \sqrt{2}, -(1/\sqrt{2}), -(1/2)\} \\ & \{(3/8) + (3/8\sqrt{2}), \sqrt{2}, (1/\sqrt{2}), (1/2)\} \end{aligned} $$
and
$$ \text{Plugging in these data sets, we get the linear system} $$

$$ \beta_0 - \sqrt{2}\beta_1 + \frac{1}{\sqrt{2}}\beta_2 - \frac{1}{2}\beta_3 = -\frac{1}{8} + \frac{1}{8\sqrt{2}} \quad (3) $$

$$ \beta_0 - \sqrt{2}\beta_1 - \frac{1}{\sqrt{2}}\beta_2 + \frac{1}{2}\beta_3 = \frac{3}{8} - \frac{3}{8\sqrt{2}} \quad (4) $$

$$ \beta_0 + \sqrt{2}\beta_1 - \frac{1}{\sqrt{2}}\beta_2 - \frac{1}{2}\beta_3 = -\frac{1}{8} - \frac{1}{8\sqrt{2}} \quad (5) $$

$$ \beta_0 + \sqrt{2}\beta_1 + \frac{1}{\sqrt{2}}\beta_2 + \frac{1}{2}\beta_3 = \frac{3}{8} + \frac{3}{8\sqrt{2}} \quad (6) $$

To estimate $\beta$ we will use the popular ‘least squares’ method, which minimises the residual sum of squares $\sum_{i=1}^N (y_i - x_i\beta_i)^2$. If $X$ is positive definite (and in turn has full rank), we can obtain a unique solution for the best fit $\hat{\beta}$, which is $(X^T X)^{-1} X^T y$. It can happen that all the columns of $X$ are not linearly independent and by extension $X$ is not full rank. This kind of situation might occur if two or more of the feature variables are perfectly correlated. Then $X^T X$ would be singular and $\hat{\beta}$ would not be uniquely defined. Nevertheless, there exist techniques like ‘filtering’ to resolve the non-unique representations by reducing the redundant features. Rank deficiencies might also occur if the number of features $p$ exceeds the number of data sets $N$. If we estimate such models using ‘regularisation’, then redundant columns should not be left out. The regulation takes care of the singularities. More importantly, the final prediction might depend on which columns are left out [22].

Equations (3)–(6) may be expressed in the matrix notation as

$$ \begin{bmatrix} -\sqrt{2} & 1 & \frac{1}{\sqrt{2}} & -\frac{1}{2} \\ -\sqrt{2} & 1 & -\frac{1}{\sqrt{2}} & \frac{1}{2} \\ -\sqrt{2} & -1 & \frac{1}{\sqrt{2}} & \frac{1}{2} \\ \sqrt{2} & 1 & \frac{1}{\sqrt{2}} & \frac{1}{2} \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_0 \\ \beta_2 \\ \beta_3 \end{bmatrix} = \begin{bmatrix} -\frac{1}{8} + \frac{1}{8\sqrt{2}} \\ \frac{3}{8} - \frac{3}{8\sqrt{2}} \\ \frac{1}{8} + \frac{1}{8\sqrt{2}} \\ \frac{3}{8} + \frac{3}{8\sqrt{2}} \end{bmatrix} \quad (7) $$

Note that unlike common convention, our representation of $X$ does not contain a column full of 1's corresponding the bias term. This representation is used simply because of the convenient form that we obtain for $X^T X$. The final result remains unaffected as long as $y = X\beta$ represents the same linear system.

Now

$$ X^T X = \frac{1}{4} \begin{bmatrix} 15 & 9 & 5 & -3 \\ 9 & 15 & 3 & -5 \\ 5 & 3 & 15 & -9 \\ -3 & -5 & -9 & 15 \end{bmatrix} \quad (8) $$

and

$$ X^T y = \begin{bmatrix} \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \\ \frac{1}{2} \end{bmatrix}. \quad (9) $$

Thus, we need to solve for $\hat{\beta}$ from $X^T X \hat{\beta} = X^T y$ [12].

## 4 Quantum circuit

Having discussed the general idea behind the HHL algorithm in Section 2 and its possible application in drastically speeding up multiple regression in Section 3, we now move on to the quantum circuit design meant to solve the 4 × 4 linear system which we encountered in Section 3, i.e. $X^T X \hat{\beta} = X^T y$. For sake of convenience, we will now denote $X^T X$ with $A$, $\hat{\beta}$ with $x$ and $X^T y$ with $b$. The circuit requires only 7 qubits, with 4 qubits in the ‘clock register’, 2 qubits in the ‘input register’ and the remaining 1 as an ‘ancilla’ qubit. At this point, it is imperative to mention that we specifically chose the form of the regression data points in the previous section such that $A$ turns out to be Hermitian, has four distinct eigenvalues of the form $\lambda_i = 2^{i-1}$ and $b$ has a convenient form which can be efficiently prepared by simply using two Hadamard gates

$$ A = X^T X = \frac{1}{4} \begin{bmatrix} 15 & 9 & 5 & -3 \\ 9 & 15 & 3 & -5 \\ 5 & 3 & 15 & -9 \\ -3 & -5 & -9 & 15 \end{bmatrix} \quad (10) $$

$A$ is a Hermitian matrix with eigenvalues $\lambda_1 = 1$, $\lambda_2 = 2$, $\lambda_3 = 4$ and $\lambda_4 = 8$. The corresponding eigenvectors encoded in quantum states $|u_j\rangle$ may be expressed as

$$ |u_1\rangle = -|00\rangle - |01\rangle - |10\rangle + |11\rangle \quad (11) $$

$$ |u_2\rangle = +|00\rangle + |01\rangle - |10\rangle + |11\rangle \quad (12) $$

$$ |u_3\rangle = +|00\rangle - |01\rangle + |10\rangle + |11\rangle \quad (13) $$

$$ |u_4\rangle = -|00\rangle + |01\rangle + |10\rangle + |11\rangle \quad (14) $$

Also

$$ b = \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{bmatrix}^T $$

can be written as $\sum_{j=1}^{j=4} \beta_j |u_j\rangle$ where each $\beta_j = (1/2)$.

We will now trace through the quantum circuit in Fig. 2. $|q_0\rangle$ and $|q_1\rangle$ are the input register qubits which are initialised to a combined quantum state

$$ |b\rangle = \frac{1}{2} |00\rangle + \frac{1}{2} |01\rangle + \frac{1}{2} |10\rangle + \frac{1}{2} |11\rangle \quad (15) $$

which is basically the state-encoded format of $b$. This is followed by the quantum phase estimation step which involves a Walsh–Hadamard transform on the clock register qubits $|j_0\rangle$, $|j_1\rangle$, $|j_2\rangle$, $|j_3\rangle$, a clock register controlled unitary gates $U^{2^0}$, $U^{2^1}$, $U^{2^2}$ and $U^{2^3}$, where $U = \exp(iAt/16)$ and an inverse quantum Fourier transform on the clock register. As discussed in Section 2, this step would produce the state $(1/2)|0001\rangle^C |u_1\rangle^I + (1/2)|0010\rangle^C |u_2\rangle^I + (1/2)|0100\rangle^C |u_3\rangle^I + (1/2)|1000\rangle^C |u_4\rangle^I$
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Quantum circuit for solving a 4 × 4 linear system **Ax = b**. Here $f = \exp(iAt/16)$ and $\tau = \pi/2^r$. The top qubit ($|s\rangle$) is the ancilla qubit. The four qubits in the middle ($|j_0\rangle$, $|j_1\rangle$, $|j_2\rangle$ and $|j_3\rangle$) stand for the clock register C. The two qubits at the bottom ($|q_0\rangle$ and $|q_1\rangle$) form the input register 1 and two Hadamard gates are applied on them to initialise the state $|b\rangle$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** GLOA is employed to approximately decompose the $\exp(iAt/2^k)$ gates in the Hamiltonian simulation step into elementary quantum gates. The decomposition is not unique. The specific gate decompositions shown in the diagram were taken by Cao et al. [17] and the angle shifts were corrected using the `scipy.optimize.minimize` module</figcaption></figure>

, which is essentially the same as $\sum_{j=1}^{N} \beta_j |u_j\rangle \langle(\tilde{\lambda}_{j\tau_0}/2\pi)\rangle^C$, assuming $\tau_0 = 2\pi$. Also, in this specific example $|\tilde{\lambda}_j\rangle = |\lambda_j\rangle$, since the 4 qubits in the clock register are sufficient to accurately and precisely represent the 4 eigenvalues in binary. As far as the endianness of the combined quantum states is concerned, we must keep in mind that in our circuit $|q_0\rangle$ is the *most significant* qubit and $|q_3\rangle$ is the *least significant* qubit.

Next is the $R(\tilde{\lambda}^{-1})$ rotation step. We make use of an ancilla qubit $|s\rangle$ (initialised in the state $|0\rangle$), which gets phase shifted depending upon the clock register's state. Let us take an example clock register state $|0100\rangle^C = |0\rangle_{q_0}^C \otimes |1\rangle_{q_1}^C \otimes |0\rangle_{q_2}^C \otimes |0\rangle_{q_3}^C$ (binary representation of the eigenvalue corresponding to $|u_3\rangle$, that is 4). In this combined state, $|q_1\rangle$ is in the state $|1\rangle$ while $|q_0\rangle$, $|q_2\rangle$ and $|q_3\rangle$ are all the in state $|0\rangle$. This state will only trigger the $R_y(2\pi/2^r)$ rotation gate, and none of the other phase shift gates. Thus, we may say that the smallest eigenvalue states in C cause the largest ancilla rotations. Using linearity arguments, it is clear that if the clock register state had instead been $|b\rangle$, as in our original example, the final state generated by this rotation step would be $\sum_{j=1}^{N} \beta_j |u_j\rangle \langle\tilde{\lambda}_j|^C ((1-C^2/\tilde{\lambda}_j^2)^{1/2} |0+(C/\tilde{\lambda}_j)\rangle)^S$ where $C = 8\pi/2^r$. For this step, an a priori knowledge of the eigenvalues of A was necessary to design the gates. For more general cases of eigenvalues, one may refer to [17].

Then, as elaborated in Section 2, the inverse phase estimation step essentially reverses the quantum phase estimation step. The state produced by this step, conditioned on obtaining $|1\rangle$ in ancilla

is $(8\pi/2^r) \sum_{j=1}^{j=4} \frac{(1/2)}{2^{j-1}} |u_j\rangle$. Upon writing in the standard basis and normalising, it becomes $(1/\sqrt{340})(-|00\rangle + 7|01\rangle + 11|10\rangle + 13|11\rangle)$. This is proportional to the exact solution of the system $\mathbf{x} = (1/32)[-1 \ 7 \ 11 \ 13]^T$.

## 5 Simulation

We simulated the quantum circuit in Fig. 2 using the Qiskit Aer QasmSimulator backend [23], which is a noisy quantum circuit simulator backend. One of the main hurdles while implementing the quantum program was dealing with the Hamiltonian simulation step, i.e. implementing the controlled unitary $\mathbf{U} = e^{iAt}$. Taking a cue from the Cao et al. paper [17] which employed the group leaders optimisation algorithm (GLOA) [24], we approximately decomposed the $\mathbf{U}$ gate into elementary quantum gates, as shown in Fig. 3. It is however important to keep in mind that using GLOA to decompose the $\mathbf{U}$ is useful only when the matrix exponential $e^{iAt}$ is readily available. Let us call the resulting approximated unitary $\tilde{\mathbf{U}}$. The parameters of the $R_x$ and $R_{zz}$ gates given in [17] were refined using the `scipy.optimize.minimize` function [25], to minimise the Hilbert-Schmidt norm of $\mathbf{U} - \tilde{\mathbf{U}}$ (which is a measure of the relative error between $\mathbf{U}$ and $\tilde{\mathbf{U}}$). The `scipy.optimize.minimize` function makes use of the quasi-Newton algorithm of Broyden, Fletcher, Goldfarb and Shanno [26] by default. Also, we noticed that it is necessary to use a controlled-Z gate instead of a single qubit Z gate (as in [17]).

A sample output of the Qiskit code has been shown in Fig. 4.
Qiskit Simulation – Sample Output:
(Unnormalized)

Predicted solution:
[-1 7 11 13]

Simulated experiment solution:
[-0.8425 6.9604 10.9980 13.0341]

Error in found solution: 0.1660

<figure><img src="image_5.png" /><figcaption>Fig. 5 Blue bars represent the predicted state amplitudes for the solution vector |x⟩. The red bars indicate the simulated experiment results using Qiskit Aer QasmSimulator (a noisy quantum simulator backend) for the state amplitudes. The states have not been normalised, so that the deviation from the solution of the original linear equation system is prominent</figcaption></figure>

All the results have been rounded to four decimal places. The 'Error in found solution' is in essence the 2-norm of the difference between the exact solution and the output solution. We have neglected normalisation and constant factors like (1/32) in the displayed solutions.

It is clear from the low error value that the gate decomposition we used for $\bar{U}$ helps to approximately replicate the ideal circuit involving $U$. Also, the difference between the predicted and simulated results (with shots = $10^5$) has been shown in Fig. 5.

One should remember that using genetic algorithms like GLOA for the Hamiltonian simulation step is a viable technique only when cost of computing the matrix exponential $e^{iAt}$ is substantially lower than the cost of solving the corresponding linear system classically. A popular classical algorithm for calculating matrix exponentials is the 'the scaling and squaring algorithm' by Al-Mohy and Higham [27–29], whose cost is $O(n^3)$, and which is generally used along with the Padé approximation by Matlab [30] and SciPy [25]. However, this algorithm is mostly used for small dense matrices. For large sparse matrices, better approaches exist. For instance, in the Krylov space approach [31, 32], an Arnoldi algorithm is used whose cost is in the ballpark of $O(mn)$, where $n$ is the matrix size and $m$ is the number of Krylov vectors which need to be computed. In general, for large sparse matrices, GLOA may be useful but that decision needs to be made on a case-by-case basis depending on the properties of the matrix $A$.

## 6 Discussion and conclusions

We have noticed that any multiple linear regression problem can be reduced to an equivalent QLSP. This allows us to utilise the general quantum speed-up techniques like the HHL algorithm. However, for HHL we need low-error low-cost circuits for the Hamiltonian

simulation and controlled rotation steps to get accurate results. There already exist well-defined deterministic approaches like the celebrated Solovay–Kitaev and Trotter–Suzuki algorithms [33–36] that can help to decompose the Hamiltonian simulation unitary $U$ up to arbitrary accuracy. However in most cases, they provide neither minimum cost nor efficiently generable gate sequences for engineering purposes. The Solovay–Kitaev algorithm, which scales as $O(\log(1/e))$, takes advantage of the fact that quantum gates are elements of the unitary group $U(d)$, which is a Lie group and a smooth differentiable manifold, in which one can do precise calculations regarding the geometry and distance between points. Similarly, the Trotter–Suzuki algorithm uses a product formula that can precisely approximate the Hamiltonian exponential, and which exploits the structure of commutation relations in the underlying Lie algebras. In practical engineering scenarios, however, one would prefer a low-cost gate sequence that approximates the unitary just well enough, rather than a high-cost and exact or almost exact decomposition.

This is where stochastic genetic algorithms like GLOA come into play. These are heuristic algorithms that do not have well-defined time complexities and error bounds. Given a particular matrix $A$, if the time taken to find a circuit approximation for the Hamiltonian evolution is (say) polynomial-time $O(n^c)$ for some constant $c$, the exponential speedup disappears. This remains an issue with GLOA because we still need to classically compute the $2^n \times 2^n$ unitary matrix of the evolution we are trying to approximate. The real advantage of GLOA shows up when dealing with large data sets, as the standard GLOA [17, 24] restricts the number of gates to a maximum of 20 (a detailed explanation of GLOA is available in the supplementary material). So, while the number of linear equations and the size of the data sets can be arbitrarily large, it will only ever use a limited number of gates for approximating the Hamiltonian simulation. This property of GLOA ensures that time taken for Hamiltonian simulation remains nearly independent of the size of the data sets as they grow large. Our circuit design technique will be useful for engineering purposes, particularly when the regression analysis is to be done on extensive amounts of data. Using GLOA, it is also often possible to incorporate other desired aspects of circuit design into the optimisation, like specific choices of gate sets. The maximum size of data set on which we can use this approach is only limited by the number of qubits available in the quantum processors.

It is noteworthy that regression analysis is widely used in various scientific and business applications, as it enables one to understand the relationship between variable parameters and make predictions about unknowns. Every day in the actual world, the handling of larger and larger data sets craves for faster and more efficient approaches to such analyses. It is natural to ask whether it might be possible to leverage quantum algorithms for this purpose. The answer seems to be a yes, although there remain several technical and engineering challenges; the dearth of sufficient experimental realisations is testimony to this. Our work on quantum multiple linear regression, which is a first of its kind, shows a practical circuit design method that we believe will pave the path for more viable and economic experimental implementations of quantum machine learning protocols, addressing problems in business analytics, machine learning, and artificial intelligence. In particular, this work opens a new door of possibilities for time- and cost-efficient experimental realisations of various quantum machine learning protocols, exploiting heuristic approaches like GLOA, which will make incrementally large amounts of data tractable.

## 7 Acknowledgments

The authors Sanchayan Dutta, Adrien Suau, Sagnik Dutta and Suvadeep Roy contributed equally to the work. The authors Bikash K. Behera and Prasanta K. Panigrahi supervised the project. Sanchayan Dutta thanks IISER Kolkata for providing hospitality and support during the period of the Summer Student Research Program (2018). Adrien Suau acknowledges the support provided by CERFACS. Suvadeep Roy and Sagnik Dutta are obliged to INSPIRE, MHRD for their support. B.K.B. acknowledges the