Randomized Row and Column Iterative Methods with a Quantum Computer

TL;DR

This paper develops quantum algorithms for classical iterative linear solvers, achieving exponential speedups over traditional methods by leveraging block-encoding and efficient state preparation.

Contribution

It introduces fast quantum implementations of the Kaczmarz and coordinate descent methods, significantly improving their computational efficiency.

Findings

Quantum algorithms achieve exponential speedup over classical methods.

Complexity is nearly linear in the number of steps.

Assumes efficient quantum state preparation for rows and columns.

Abstract

We consider the quantum implementations of the two classical iterative solvers for a system of linear equations, including the Kaczmarz method which uses a row of coefficient matrix in each iteration step, and the coordinate descent method which utilizes a column instead. These two methods are widely applied in big data science due to their very simple iteration schemes. In this paper we use the block-encoding technique and propose fast quantum implementations for these two approaches, under the assumption that the quantum states of each row or each column can be efficiently prepared. The quantum algorithms achieve exponential speed up at the problem size over the classical versions, meanwhile their complexity is nearly linear at the number of steps.

Equations

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\left(\frac{{\mathbf{a}}_{i_k}^T{\mathbf{x}}_k}{\Vert {\mathbf{a}}_{i_k}\Vert }-\frac{b_{i_k}}{\Vert {\mathbf{a}}_{i_k}\Vert }\right)\frac{{\mathbf{a}}_{i_k}}{\Vert {\mathbf{a}}_{i_k}\Vert },$$

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{c}}_{j_k}^T(\mathbf{b}-A{\mathbf{x}}_k)}{\Vert {\mathbf{c}}_{j_k}{\Vert }^2}{\mathbf{e}}_{j_k},$$

$$\begin{array}[]{lll}\vspace{.2cm}{\bf x}_{k+1}&=&\displaystyle{\bf x}_{k}+\frac{{\bf c}_{j_{k}}^{T}{\bf r}_{k}}{\|{\bf c}_{j_{k}}\|^{2}}{\bf e}_{j_{k}},\\ {\bf r}_{k+1}&=&\displaystyle{\bf r}_{k}-\frac{{\bf c}_{j_{k}}{\bf c}_{j_{k}}^{T}}{\|{\bf c}_{j_{k}}\|^{2}}{\bf r}_{k}.\end{array}$$

$$|{\bf x}_{k+1}\rangle\varpropto\|{\bf x}_{k}\|\Big{(}I_{n}-|{\bf a}_{i_{k}}\rangle\langle{\bf a}_{i_{k}}|\Big{)}|{\bf x}_{k}\rangle+b_{i_{k}}|{\bf a}_{i_{k}}\rangle.$$

$$U_{t}=\left[\begin{array}[]{cc}\vspace{.2cm}I_{n}-|{\bf a}_{t}\rangle\langle{\bf a}_{t}|&~{}~{}|{\bf a}_{t}\rangle\langle{\bf a}_{t}|\\ |{\bf a}_{t}\rangle\langle{\bf a}_{t}|&~{}~{}I_{n}-|{\bf a}_{t}\rangle\langle{\bf a}_{t}|\\ \end{array}\right]=I_{2}\otimes(I_{n}-|{\bf a}_{t}\rangle\langle{\bf a}_{t}|)+X\otimes|{\bf a}_{t}\rangle\langle{\bf a}_{t}|,$$

$$U_t=(I_2\otimes V_t)(I_2\otimes (I_n-|0⟩⟨0|)+X\otimes |0⟩⟨0|)(I_2\otimes V_t^{†}).$$

$$|X⟩=\sqrt{p}|0⟩|{\mathbf{x}}_k⟩+\sqrt{1-p}|1⟩|\cdots ⟩.$$

$$|\psi\rangle={\rm SWAP}_{1,2}\Big{(}\beta|0\rangle|X\rangle+\gamma|1\rangle|0\rangle|{\bf a}_{t}\rangle\Big{)}=|0\rangle\Big{(}\beta\sqrt{p}\,|0\rangle|{\bf x}_{k}\rangle+\gamma|1\rangle|{\bf a}_{t}\rangle\Big{)}+\beta\sqrt{1-p}\,|1\rangle|0\rangle|\cdots\rangle.$$

$$U_{t}\Big{(}\beta\sqrt{p}\,|0\rangle|{\bf x}_{k}\rangle+\gamma|1\rangle|{\bf a}_{t}\rangle\Big{)}=|0\rangle\otimes\Big{(}\beta\sqrt{p}(I_{n}-|{\bf a}_{t}\rangle\langle{\bf a}_{t}|)|{\bf x}_{k}\rangle+\gamma|{\bf a}_{t}\rangle\Big{)}+\beta\sqrt{p}\langle{\bf a}_{t}|{\bf x}_{k}\rangle|1\rangle|{\bf a}_{t}\rangle.$$

$$|X_k⟩=\frac{\Vert {\mathbf{x}}_k\Vert }{{\mu }_k}|0{⟩}^{\otimes k}\otimes |{\mathbf{x}}_k⟩+|0^{\perp }{⟩}^{\otimes k}|\cdots ⟩.$$

$$|X_{k+1}⟩=\frac{\Vert {\mathbf{x}}_{k+1}\Vert }{{\mu }_{k+1}}|0{⟩}^{\otimes (k+1)}\otimes |{\mathbf{x}}_{k+1}⟩+|0^{\perp }{⟩}^{\otimes (k+1)}|\cdots ⟩.$$

$$\begin{array}[]{lll}\vspace{.2cm}&&\displaystyle I_{2}^{\otimes k}\otimes U_{t_{k}}\Big{(}\frac{\beta_{t_{k}}\|{\bf x}_{k}\|}{\mu_{k}}|0\rangle^{\otimes k}|0\rangle|{\bf x}_{k}\rangle+\gamma_{t_{k}}|0\rangle^{\otimes k}|1\rangle|{\bf a}_{t_{k}}\rangle+|0^{\bot}\rangle^{\otimes k}|0\rangle|\cdots\rangle\Big{)}\\ \vspace{.2cm}&=&\displaystyle|0\rangle^{\otimes(k+1)}\otimes\Big{(}\frac{\beta_{t_{k}}\|{\bf x}_{k}\|}{\mu_{k}}\Big{(}I_{n}-|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|\Big{)}|{\bf x}_{k}\rangle+\gamma_{t_{k}}|{\bf a}_{t_{k}}\rangle\Big{)}+|0^{\bot}\rangle^{\otimes(k+1)}|\cdots\rangle\\ &=&\displaystyle\frac{\beta_{t_{k}}}{\mu_{k}}|0\rangle^{\otimes(k+1)}\otimes\Big{(}\|{\bf x}_{k}\|\Big{(}I_{n}-|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|\Big{)}|{\bf x}_{k}\rangle+b_{t_{k}}|{\bf a}_{t_{k}}\rangle\Big{)}+|0^{\bot}\rangle^{\otimes(k+1)}|\cdots\rangle,\end{array}$$

$${\mu }_k^2=1+\sum _{i=1}^{k-1}{b}_{t_i}^2.$$

$${\mu }_{k+1}^2=\frac{{\mu }_k^2}{{\beta }_{t_k}^2}={\mu }_k^2+b_{t_k}^2.$$

$${\mu }_k^2=1+\sum _{i=1}^{k-1}{b}_{t_i}^2.$$

$$\mathbb{E}[{\mu }_k^2]=1+(k-1)\mathbb{E}[b_{t_i}^2]=1+\frac{k-1}{n}\Vert \mathbf{b}{\Vert }_2^2.$$

$$\begin{array}[]{rll}\vspace{.2cm}{\bf x}_{k+1}&=&\displaystyle{\bf x}_{k}+\frac{{\bf a}_{t_{k}}^{T}{\bf r}_{k}}{\|{\bf a}_{t_{k}}\|^{2}}{\bf e}_{t_{k}},\\ {\bf r}_{k+1}&=&\displaystyle\Big{(}I_{n}-\frac{{\bf a}_{t_{k}}{\bf a}_{t_{k}}^{T}}{\|{\bf a}_{t_{k}}\|^{2}}\Big{)}{\bf r}_{k}.\end{array}$$

$$\begin{array}[]{rll}\vspace{.2cm}|{\bf x}_{k+1}\rangle&\varpropto&\displaystyle\|{\bf x}_{k}\|\,|{\bf x}_{k}\rangle+\|{\bf r}_{k}\|\,|{t_{k}}\rangle\langle{\bf a}_{t_{k}}|{\bf r}_{k}\rangle,\\ |{\bf r}_{k+1}\rangle&\varpropto&\displaystyle\|{\bf r}_{k}\|\,(I_{n}-|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|)|{\bf r}_{k}\rangle.\end{array}$$

$$|R_k⟩=\Vert {\mathbf{r}}_k\Vert |0{⟩}^{\otimes k}|{\mathbf{r}}_k⟩+|0^{\perp }{⟩}^{\otimes k}|\cdots ⟩.$$

$$\begin{array}[]{lll}\vspace{.2cm}(I_{2}^{\otimes k}\otimes U_{t_{k}}){\rm SWAP}_{1,k+1}|0\rangle|R_{k}\rangle&=&|0\rangle^{\otimes(k+1)}\|{\bf r}_{k}\|(I_{n}-|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|)|{\bf r}_{k}\rangle+|0^{\bot}\rangle^{\otimes(k+1)}|\cdots\rangle\\ \vspace{.2cm}&=&\|{\bf r}_{k+1}\|\,|0\rangle^{\otimes(k+1)}|{\bf r}_{k+1}\rangle+|0^{\bot}\rangle^{\otimes(k+1)}|\cdots\rangle\\ &=&|R_{k+1}\rangle.\end{array}$$

$$|{\mathbf{x}}_{k+1}⟩\propto \Vert {\mathbf{x}}_k\Vert |{\mathbf{x}}_k⟩+\Vert {\mathbf{r}}_k\Vert |t_k⟩⟨t_k|S_{t_k}|{\mathbf{r}}_k⟩.$$

$$|{\tilde{X}}_k⟩=\frac{\Vert {\mathbf{x}}_k\Vert }{\mu }|0⟩|{\mathbf{x}}_k⟩+|0^{\perp }⟩|\cdots ⟩,$$

$$|{\tilde{R}}_k⟩=\Vert {\mathbf{r}}_k\Vert |0⟩S_{t_k}|{\mathbf{r}}_k⟩+|0^{\perp }⟩|\cdots ⟩.$$

$$|{\varphi }_1⟩=\beta |00⟩|{\tilde{X}}_k⟩+\gamma |10⟩|{\tilde{R}}_k⟩,$$

$$W_{t}=\left[\begin{array}[]{ccc}\vspace{.2cm}I_{n}&0&0\\ \vspace{.2cm}0&I_{n}-|t\rangle\langle t|&|t\rangle\langle t|\\ 0&|t\rangle\langle t|&I_{n}-|t\rangle\langle t|\\ \end{array}\right].$$

$$\begin{array}[]{lll}\vspace{.2cm}|\phi_{2}\rangle&=&{\rm SWAP}_{2,3}W_{t_{k}}|\phi_{1}\rangle\\ \vspace{.2cm}&=&{\rm SWAP}_{2,3}(\beta|00\rangle|\widetilde{X}_{k}\rangle+\gamma|01\rangle|t_{k}\rangle\langle t_{k}|\widetilde{R}_{k}\rangle+|10\rangle|\cdots\rangle)\\ &=&\displaystyle|00\rangle\Big{(}\frac{\beta\|{\bf x}_{k}\|}{\mu}|0\rangle|{\bf x}_{k}\rangle+\gamma\|{\bf r}_{k}\|\,|1\rangle|t_{k}\rangle\langle t_{k}|S_{t_{k}}|{\bf r}_{k}\rangle\Big{)}+|0^{\bot}\rangle^{\otimes 2}|\cdots\rangle.\end{array}$$

$$|\phi_{3}\rangle=(I_{4}\otimes G_{k}\otimes I_{n})|\phi_{2}\rangle=|000\rangle\Big{(}\frac{c\beta\|{\bf x}_{k}\|}{\mu}|{\bf x}_{k}\rangle+s\gamma\|{\bf r}_{k}\|\,|t_{k}\rangle\langle t_{k}|S_{t_{k}}|{\bf r}_{k}\rangle\Big{)}+|0^{\bot}\rangle^{\otimes 3}|\cdots\rangle.$$

$$\Vert {\mathbf{x}}_k\Vert |{\mathbf{x}}_k⟩+\Vert {\mathbf{r}}_k\Vert |t_k⟩⟨t_k|S_{t_k}|{\mathbf{r}}_k⟩=\Vert {\mathbf{x}}_{k+1}\Vert |{\mathbf{x}}_{k+1}⟩.$$

$$G_{k}=\frac{1}{\sqrt{k+2}}\left[\begin{array}[]{cc}\vspace{.2cm}\sqrt{k+1}&1\\ -1&\sqrt{k+1}\\ \end{array}\right].$$

$$\begin{array}[]{lll}\vspace{.2cm}|X_{k}\rangle&=&\displaystyle\frac{\|{\bf x}_{k}\|}{k+1}\,|0\rangle^{\otimes 2k}\otimes|{\bf x}_{k}\rangle+|0^{\bot}\rangle^{\otimes 2k}|\cdots\rangle,\\ |R_{k}\rangle&=&\|{\bf r}_{k}\|\,|0\rangle^{\otimes k}|{\bf r}_{k}\rangle+|0^{\bot}\rangle^{\otimes k}|\cdots\rangle.\end{array}$$

Full text

Keywords. Kaczmarz method, coordinate descent method, block-encoding, randomized algorithms, quantum iterative algorithms

1 Introduction

The classical solvers for a linear system of equations $A{\bf x}={\bf b}$ are generally categorized into two types: direct methods and iterative methods. The latter is usually more practical in the realm of large-scale system of equations. Among the iterative methods, the Kaczmarz method and the coordinate descent method are popular due to the simplicity and efficiency. The Kaczmarz method was first discovered in 1937 [10], and was rediscovered in the field of image reconstruction by Gordon, Bender and Herman in 1970 [8] under the appellation algebraic reconstruction technique (ART). The Kaczmarz method uses a row of coefficient matrix in each iteration, while the coordinate descent method utilizes a column instead. These two methods seek to solve different problems: the coordinate descent method converges to a least squares solution generally; the Kaczmarz method calculates a minimum-norm solution for a consistent system of equations and exhibits cyclic convergence for an inconsistent problem [5]. They can be generalized to many variants (see [16] and the references therein). The advantage of these two methods lies in the fact that at each iteration they only need access to an individual row (or column) rather than the entire coefficient matrix. Due to the simplicity, they have numerous applications in the fields ranging from computer tomography to image processing and digital signal processing, especially the big data science.

We review the iteration schemes in the following. Assume that $A$ is an $n$-by-$n$ matrix. For any $1\leq i\leq n$, denote the $i$-th row of $A$ as ${\bf a}_{i}^{T}$, and the $i$-th component of ${\bf b}$ as $b_{i}$. Let ${\bf x}_{0}$ be an arbitrary initial approximation to the solution of $A{\bf x}={\bf b}$. For $k\geq 0$, randomly choose an $i_{k}\in\{1,\ldots,n\}$, the Kaczmarz iteration updates the solution ${\bf x}_{k}$ as follows:

[TABLE]

which is equivalent to the Gauss-Seidel method on $AA^{T}{\bf y}={\bf b}$, where $A^{T}{\bf y}={\bf x}$. Geometrically, ${\bf x}_{k+1}$ is the orthogonal projection of ${\bf x}_{k}$ onto the hyperplane ${\bf a}_{i_{k}}^{T}{\bf x}=b_{i_{k}}$. In each iteration step, only one row of the coefficient matrix is needed. It is also called the row-action method. Another names, such as component-solution method, cyclic projection or successive projection, are used on certain occasions.

Correspondingly, we have the column-action method. Let ${\bf c}_{j}$ be the $j$-th column of $A$, and ${\bf e}_{j}$ the $j$-th column of the unit matrix. The column-action method reads

[TABLE]

where $j_{k}$ is a random number from $\{1,\ldots,n\}$. It is equivalent to the randomized Gauss-Seidel method on $A^{T}A{\bf x}=A^{T}{\bf b}$. In each iteration step the approximate solution is obtained by a minimization in one coordinate direction, then it is also called the coordinate descent method. Define ${\bf r}_{k}={\bf b}-A{\bf x}_{k}$. This iteration method can be reexpressed as

[TABLE]

Suppose that the size of the problem is $n$ and the number of required iterations to solve the problem is $T$, then the complexity of classical iteration algorithm is usually polynomial at $n$ and linear at $T$. Due to the quantum no-clone theorem, an iteration algorithm is usually not easy to implement in a quantum computer.

For instance, in [13], Rebentrost et al. proposed a quantum gradient and Newton’s method to solve polynomial optimization. Compared to the classical gradient or Newton’s method, this quantum algorithm achieves exponential speedup at $n$. However, the complexity exponentially depends on $T$. The main reason is as follows, taking Newton’s method as an example. One critical step of classical Newton’s method is to solve a linear system. The coefficient matrix (i.e., the Hessian matrix) $H$ of the linear system depends on ${\bf x}_{k}$, the computed result in step $k$. In a quantum computer, we only have the quantum state of ${\bf x}_{k}$, which is a unknown quantum state. To solve the linear system by HHL algorithm, one critical step is the Hamiltonian simulation of $H$. In [13], a similar idea to quantum principal component analysis [12] is used to compute the Hamiltonian simulation of $H$ by viewing it as a unknown density matrix. Together with quantum phase estimation, the linear system can be solved efficiently. However, the Hamiltonian simulation of a unknown density operator requires $O(t^{2}/\epsilon)$ copies of this density operator, where $t$ is the evolution time and $\epsilon$ is the precision. Since $H$ depends on $|{\bf x}_{k}\rangle$, we need to prepare at least $O(t^{2}/\epsilon)$ copies of $|{\bf x}_{k}\rangle$. In other words, with at least $O(t^{2}/\epsilon)$ copies of $|{\bf x}_{k}\rangle$, we can prepare $|{\bf x}_{k+1}\rangle$. As a result, we need exponentially copies of the initial state to prepare $|{\bf x}_{T}\rangle$.

In [11], Kerenidis and Prakash considered a quantum version of interior-point method to solve semi-definite programming and linear programming. One critical step of classical interior-point method is also to solve a linear system. Different from the idea used above, they output the classical information of $|{\bf x}_{k}\rangle$ at each step of iteration. By doing so $|{\bf x}_{k}\rangle$ becomes a known quantum state. To read out the classical information of a quantum state by quantum tomography requires at least $O(n)$ steps. Therefore, the complexity of their quantum algorithm is polynomial at $n$. The dependence on $T$ is polynomial due to the error caused in quantum tomography in each step.

These are two typical examples of quantum iteration algorithms. They cannot outperform classical iteration algorithms both at $n$ and $T$. The quantum iteration algorithm is a simulation of classical iteration algorithm in a quantum computer, so it seems especially hard to achieve a speedup at $T$. However, we still expect the quantum iteration algorithm to achieve high speedup at $n$ and has a reasonable dependence on $T$ meanwhile. Therefore, an ideal simulation of iteration algorithm in a quantum computer should have the complexity polynomial-logarithm at $n$ and linear or polynomial at $T$. In this paper, we will give such implementations to the Kaczmarz method and the coordinate descent method. The main idea is to use block-encoding technique [2] to implement the procedures (1), (3) by unitary operators. These unitary operators are explicitly constructed and efficiently implemented in the quantum computer. Moreover, they can help us overcome the calculations of the inner product of quantum states. As a result, we can show that the Kaczmarz method and the coordinate descent method can be implemented in a quantum computer in time $O(T(\log n))$.

Throughout this paper, we use the following notations.

1. i.

A quantum state of the form $\sum_{i=0}^{2^{k}-1}\alpha_{i}|i\rangle|\psi_{i}\rangle$ will be simply denoted by $\alpha_{0}|0\rangle^{\otimes k}|\psi_{0}\rangle+|0^{\bot}\rangle^{\otimes k}|\cdots\rangle$, when we only concern about $|\psi_{0}\rangle$ and neglect the garbage state. 2. ii.

We denote ${\rm SWAP}_{i,j}$ as the swap operator that swaps the $i$-th qubit and the $j$-th qubit.

2 The quantum Kaczmarz algorithm

Assume that the quantum state of ${\bf a}_{t}$ can be prepared efficiently in the quantum computer, such as by qRAM [7]. So there is an efficiently implemented unitary operator $V_{t}$ such that $V_{t}|0\rangle=|{\bf a}_{t}\rangle$. Based on the iteration formula (1), without loss of generality, we suppose that $\|{\bf a}_{t}\|=1$ for all $t$. By equation (1), we have

[TABLE]

For any row index $t$, define a unitary operator

[TABLE]

where $X$ is the Pauli-X matrix. That is,

[TABLE]

By qRAM assumption, $V_{t}$ is efficiently implemented in the quantum computer, so is $U_{t}$.

By equation (5), $U_{t}$ can be viewed as a control operator: if the second register is $|{\bf a}_{t}\rangle$, then apply $X$ to the first register; if the second register lies in the orthogonal complement space of $|{\bf a}_{t}\rangle$, then do nothing to the first register. The general architecture of these kind of unitaries was studied in [6].

The basic idea of the quantum implementation of Kaczmarz iteration is as follows: Suppose that we have the following state that contains the quantum information of $|{\bf x}_{k}\rangle$

[TABLE]

Let $\beta^{2}+\gamma^{2}=1$, then we can prepare

[TABLE]

As to the first term, direct calculation yields that

[TABLE]

Properly choosing the parameters $\beta,\gamma$, for example, $\beta=\|{\bf x}_{k}\|\delta,\gamma=b_{t}\sqrt{p}\delta$ for some $\delta$ to ensure $\beta^{2}+\gamma^{2}=1$, then the first term of $|\psi\rangle$ is a state proportional to the right hand side of equation (4).

The explicit procedure to implement Kaczmarz iteration is stated as follows.

Algorithm 1 The quantum Kaczmarz method

1:  Randomly choose a unit vector ${\bf x}_{0}$ such that its quantum state can be prepared in time $O(\log n)$. Set $k=0$ and $\mu_{k}=1$. The state can be expressed in the following general form

[TABLE]

2:  Randomly choose a $t_{k}\in\{1,\ldots,n\}$. Define $\beta_{t_{k}}^{2}=\displaystyle\frac{\mu_{k}^{2}}{\mu_{k}^{2}+b_{t_{k}}^{2}}$, $\gamma_{t_{k}}^{2}=1-\beta_{t_{k}}^{2}$ and $\mu_{k+1}=\displaystyle\frac{\mu_{k}}{\beta_{t_{k}}}$.

3:  Apply $(I_{2}^{\otimes k}\otimes U_{t_{k}}){\rm SWAP}_{1,k+1}$ to $\beta_{t_{k}}|0\rangle|X_{k}\rangle+\gamma_{t_{k}}|1\rangle|0\rangle^{\otimes k}|{\bf a}_{t_{k}}\rangle$, then we obtain

[TABLE]

4:  Set $k=k+1$, and go to step 2 until converges.

In step 3, we calculate that

[TABLE]

where we use the fact that $\gamma_{t_{k}}=\sqrt{1-\beta_{t_{k}}^{2}}=\beta_{t_{k}}b_{t_{k}}/\mu_{k}$ in the last step.

Similar to the classical Kaczmarz method, algorithm 2 is also simple to implement in a quantum computer. Let ${\bf x}_{k}$ be the result obtained by the classical Kaczmarz method in the $k$-th step. Then the first term of $|X_{k}\rangle$ defined in equation (10) contains all the information of ${\bf x}_{k}$ precisely, i.e., $\|{\bf x}_{k}\||{\bf x}_{k}\rangle$.

Theorem 1

Assume that $|{\bf a}_{t}\rangle$ is prepared in $O(\log n)$ for any $t$. In algorithm 2, for any $k\geq 1$, we have

[TABLE]

The complexity to prepare $|X_{k}\rangle$ is $O(k\log n)$.

Proof. By definition in step 2 of algorithm 2,

[TABLE]

Since $\mu_{0}=1$, we have

[TABLE]

Assume that the complexity to prepare $|X_{k}\rangle$ is $\tau_{k}$, then the complexity for $|X_{k+1}\rangle$ in step 3 is $\tau_{k}+O(\log n)$, since the time for preparing $|a_{t_{k}}\rangle$ is $O(\log n)$. Thus, $\tau_{k+1}=\tau_{k}+O(\log n)$. Since $\tau_{0}=O(\log n)$, we have $\tau_{k}=O(k\log n)$. $\hfill{\blacksquare}$

Strohmer et al. [14] sample a row in a random fashion with probability proportional to the 2-norm of that row at each iteration, and prove an exponential expected convergence rate of randomized Kaczmarz method. For simplicity, we assume a uniform sampling in step 2. Then the expectation reads

[TABLE]

If $\|{\bf b}\|_{\infty}=O(1)$, then $\mathbb{E}[\mu_{k}]=O(\sqrt{k})$.

The classical information of the solution ${\bf x}_{k}$ is stored in the first term of $|X_{k}\rangle$. For some problems in machine learning, such as data classification, the final output extracts certain global information, rather than each component of a state. For example, to estimate the inner product between ${\bf x}_{k}$ and the vector ${\bf c}$, we can first prepare the quantum state $|{\bf c}\rangle$ of ${\bf c}$, then apply the swap test to estimate the inner product between $|X_{k}\rangle$ and $|0\rangle^{\otimes k}|{\bf c}\rangle$. This returns an $\epsilon^{\prime}$-approximate of $\|{\bf x}_{k}\|\langle{\bf x}_{k}|{\bf c}\rangle/\mu_{k}$ in time $O(k(\log n)/\epsilon^{\prime})$. Thus, by setting $\epsilon^{\prime}=\epsilon/\mu_{k}\|{\bf c}\|$, we will obtain an $\epsilon$-approximate of ${\bf x}_{k}\cdot{\bf c}$ in time $O(k(\log n)\mu_{k}\|{\bf c}\|/\epsilon)$.

3 The quantum coordinate descent algorithm

With a bit abuse of notations, in the following we use ${\bf a}_{j}$ to denote the $j$-th column of $A$. Assume that the quantum state of ${\bf a}_{j}$ can be efficiently prepared. That is, there exist unitary operators $S_{j}$ such that $S_{j}^{\dagger}|j\rangle=|{\bf a}_{j}\rangle$ for any $j$.

The coordinate descent method can be implemented as follows: (1) Randomly choose an initial guess ${\bf x}_{0}$ and set ${\bf r}_{0}={\bf b}-A{\bf x}_{0}$. (2) Randomly choose a $t_{k}\in\{1,\ldots,n\}$ and update

[TABLE]

For convenience, we assume that $\|{\bf a}_{t}\|=1$ for any $t$. Using quantum states, we can rewrite (13) as

[TABLE]

Before implementing the procedure (14) in the quantum computer, we state some ideas below. Firstly, we consider the update of the residual. The basic idea is the same as algorithm 2. Suppose that the residual of the $k$-th step is encoded in the state

[TABLE]

Apply $(I_{2}^{\otimes k}\otimes U_{t_{k}}){\rm SWAP}_{1,k+1}$ to $|0\rangle|R_{k}\rangle$, then

[TABLE]

This is in fact the algorithm 2 with initial vector ${\bf r}_{0}$ and $\mu_{k}=1$ for all $k$.

Secondly, we consider the update of the approximate solution. Since $\langle t|S_{t}=\langle{\bf a}_{t}|$, we have

[TABLE]

Algorithm 2 is not applicable to the above procedure directly. Some modifications are required. The following are some basic ideas.

Suppose that the approximate solution ${\bf x}_{k}$ and the corresponding residual ${\bf r}_{k}$ are encoded in the following states, respectively

[TABLE]

and

[TABLE]

We then show how to combine them to generate a state that contains $|{\bf x}_{k+1}\rangle$.

Introduce two auxilla qubits, and prepare

[TABLE]

where $\beta^{2}+\gamma^{2}=1$. For any $t$, define

[TABLE]

Apply ${\rm SWAP}_{2,3}W_{t_{k}}$ to $|\phi_{1}\rangle$ to prepare

[TABLE]

Define G_{k}=\left[\begin{array}[]{rr}\vspace{.2cm}c&s\\ -s&c\\ \end{array}\right], where $c^{2}+s^{2}=1$. Apply $I_{4}\otimes G_{k}\otimes I_{n}$ to $|\phi_{2}\rangle$, then we obtain

[TABLE]

We can properly choose the parameters $c,s,\beta,\gamma$, such that the first term of $|\phi_{3}\rangle$ is proportional to

[TABLE]

With the preparations above, we can present the quantum coordinate descend algorithm as follows, where

[TABLE]

Algorithm 2 The quantum coordinate descend method

1:  Randomly choose a unit vector ${\bf x}_{0}$ such that its quantum state can be prepared in time $O(\log n)$. Assume that ${\bf r}_{0}={\bf b}-A{\bf x}_{0}$ has unit norm and its quantum state is prepared in time $O(\log n)$. Set $k=0$. Denote

[TABLE]

2:  Randomly choose a $t_{k}\in\{1,\ldots,n\}$.

3:  Apply $(I_{2}^{\otimes(2k+1)}\otimes G_{k}\otimes I_{n})(I_{2}^{\otimes 2k}\otimes W_{t_{k}}){\rm SWAP}_{2,2k+2}{\rm SWAP}_{1,2k+1}$ to

[TABLE] then we obtain $|X_{k+1}\rangle$.

4:  Apply $(I_{2}^{\otimes k}\otimes U_{t_{k}}){\rm SWAP}_{1,k+1}$ to $|0\rangle|R_{k}\rangle$ to generate $|R_{k+1}\rangle$.

5:  Set $k=k+1$, and go to step 2 until converges.

We explain the update of approximate solution in details. The state (27) in step 3 is denoted as $|\psi_{0}\rangle$, that is,

[TABLE]

Apply ${\rm SWAP}_{2,2k+2}{\rm SWAP}_{1,2k+1}$ to $|\psi_{0}\rangle$, we obtain

[TABLE]

Apply $I_{2}^{\otimes 2k}\otimes W_{t_{k}}$ to $|\psi_{1}\rangle$ to get

[TABLE]

Apply $I_{2}^{\otimes(2k+1)}\otimes G_{k}\otimes I_{n}$ to $|\psi_{2}\rangle$, then we get

[TABLE]

Similar to the proof of theorem 1, we obtain the following result

Theorem 2

In algorithm 3 the complexity to prepare $|X_{k}\rangle$ is $O(k\log n)$.

Remark 1

Since $\|{\bf x}_{0}\|=1$, now assume that $\|{\bf x}_{k}\|\leq k+1$ and $\|{\bf r}_{0}\|=1$, then by equation (14),

[TABLE]

Therefore, by induction the equation (26) is well-defined.

Remark 2

By definition, ${\bf r}_{0}={\bf b}-A{\bf x}_{0}$, with a suitable choice of ${\bf x}_{0}$ we can make sure that it has unit norm. Even if it does not has unit norm, algorithm 3 still works. Let $\rho$ be a parameter such that $\rho\|{\bf r}_{0}\|\leq 1$. In algorithm 3, it suffices to change $|R_{0}\rangle$ into $|\widehat{R}_{0}\rangle=\rho\|{\bf r}_{0}\|\,|0\rangle|{\bf r}_{0}\rangle+|0^{\bot}\rangle|\cdots\rangle$. By equation (14), $\|{\bf r}_{k+1}\|\,|{\bf r}_{k+1}\rangle=\|{\bf r}_{k}\|\,(I_{n}-|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|)|{\bf r}_{k}\rangle.$ Therefore, $|R_{k}\rangle$ used in algorithm 3 is simply changed into $|\widehat{R}_{k}\rangle=\rho\|{\bf r}_{k}\|\,|0\rangle^{\otimes(k+1)}|{\bf r}_{k}\rangle+|0^{\bot}\rangle^{\otimes(k+1)}|\cdots\rangle.$ Here the reason to use $k+1$ ancilla qubits is the extra one ancilla qubit introduced in $|\widehat{R}_{0}\rangle$. Since $\rho\|{\bf r}_{0}\|\leq 1$, we have $\rho\|{\bf r}_{k}\|\leq 1$ for all $k$.

Remark 3

For the update of residual we can use the linear combinations of unitaries (LCU) [3]. It contains many applications in quantum computing, such as quantum simulation [1, 3], quantum linear solver [4]. Let $U_{0},\ldots,U_{m-1}$ be $m$ unitary operators, and $\alpha_{0},\ldots,\alpha_{m-1}$ be $m$ positive real numbers. Set $s=\sum_{j}\alpha_{j}$. Assume that $V$ is a unitary operator that maps $|0\rangle^{\otimes\log m}$ to $\frac{1}{\sqrt{s}}\sum_{j}\sqrt{\alpha_{j}}|j\rangle$. Given a quantum state $|\psi\rangle$, the technique of LCU can compute $\sum_{j}\alpha_{j}U_{j}|\psi\rangle$. In our case, if we choose $U_{0}=I_{n}$, $U_{1}=I_{n}-2|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|$ and $\alpha_{0}=\alpha_{1}=1/2$, then we can prepare $(I_{n}-|{\bf a}_{t_{k}}\rangle\langle{\bf a}_{t_{k}}|)|\psi\rangle$.

4 Conclusions

The quantum implementation of a general iterative method is usually very challenging. The typical quantum iteration methods in [11, 13] cannot outperform classical iteration methods in both the problem size $n$ and the iteration number $T$. We therefore switch to two special iterative methods: the Kaczmarz method and the coordinate descent method. For solving linear systems, these two methods may not be absolutely superior to other iterative algorithms. But due to their simple structures that only one row (or column) is accessed at each iteration, they are very popular in certain areas, such as the big data science. In this paper, we show that the quantum versions of Kaczmarz method and coordinate descent method also have simple implementation structures in a quantum computer. Moreover, the efficiency is exponentially better in the problem size than their classical counterparts.

Our quantum iterative linear solvers are different from other quantum linear solvers [4, 9, 15]. The methods in this paper are independent of Hamiltonian simulation and quantum phase estimation, and also shed some lights on the quantum implementation of iterative methods. The only assumption on our algorithms is the requirement of a qRAM, by which we can efficiently extract the row or column of the coefficient matrix. One drawback of our quantum iterative linear solvers is the number of ancilla qubits. Since the iteration is generally not a unitary procedure, we need to introduce ancilla qubits to change it into a unitary one. It remains an open problem to find better ways to reduce the number of ancilla qubits.