Quantum algorithm and quantum circuit for A-Optimal Projection: dimensionality reduction

TL;DR

This paper introduces a quantum algorithm and circuit for A-Optimal Projection, achieving exponential speedup over classical methods in dimensionality reduction tasks, with efficient complexity and broad applicability.

Contribution

The paper presents a novel quantum algorithm and circuit for A-Optimal Projection, improving efficiency and removing quantum-classical interactions compared to prior classical approaches.

Findings

Exponential speedup in feature space dimensions

Quantum circuit complexity scales logarithmically

High fidelity with reduced quantum-classical interaction

Abstract

Learning low dimensional representation is a crucial issue for many machine learning tasks such as pattern recognition and image retrieval. In this article, we present a quantum algorithm and a quantum circuit to efficiently perform A-Optimal Projection for dimensionality reduction. Compared with the best-know classical algorithms, the quantum A-Optimal Projection (QAOP) algorithm shows an exponential speedup in both the original feature space dimension $n$ and the reduced feature space dimension $k$. We show that the space and time complexity of the QAOP circuit are $O\left[ {{{\log }_2}\left( {nk} /{\epsilon} \right)} \right]$ and $O[ {\log_2(nk)} {poly}\left({{\log }_2}\epsilon^{-1} \right)]$ respectively, with fidelity at least $1-\epsilon$. Firstly, a reformation of the original QAOP algorithm is proposed to help omit the quantum-classical interactions during the QAOP algorithm. Then the quantum algorithm and quantum circuit with performance guarantees are proposed. Specifically, the quantum circuit modules for preparing the initial quantum state and implementing the controlled rotation can be also used for other quantum machine learning algorithms.

Figures

Equations

$${S_{ij}}=\left\{\begin{array}[]{l}1,\quad if\quad{{\bf{x}}_{i}}\in{N_{k}}\left({{{\bf{x}}_{j}}}\right)\quad or\quad{{\bf{x}}_{j}}\in{N_{k}}\left({{{\bf{x}}_{i}}}\right)\\ 0,\quad otherwise\end{array}\right.$$

$$z={\beta }^T\mathbf{y}+{ϵ}_0$$

$$\underset{\mathbf{A}}{\min }Tr\left({\left({\mathbf{A}}^T\mathbf{X}\left(\mathbf{I}+{\lambda }_1\mathbf{L}\right){\mathbf{X}}^T\mathbf{A}+{\lambda }_2\mathbf{I}\right)}^{-1}\right)$$

$$\underset{\mathbf{A},\mathbf{B}}{\min }{\Vert \mathbf{I}-{\mathbf{A}}^T\tilde{\mathbf{X}}\mathbf{B}\Vert }^2+{\lambda }_2{\Vert \mathbf{B}\Vert }^2$$

$$\frac{\partial \varphi }{\partial {\mathbf{B}}^T}=0\Rightarrow \mathbf{B}={\left({\tilde{\mathbf{X}}}^T\mathbf{A}{\mathbf{A}}^T\tilde{\mathbf{X}}+{\lambda }_2\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^T\mathbf{A}$$

$$\frac{\partial \varphi }{\partial \mathbf{A}}=0\Rightarrow \mathbf{A}={\left(\tilde{\mathbf{X}}\mathbf{B}{\mathbf{B}}^T{\tilde{\mathbf{X}}}^T\right)}^{-1}\tilde{\mathbf{X}}\mathbf{B}$$

$${\mathbf{A}}^{\left(0\right)}=pca\left(\tilde{\mathbf{X}}\right)=\sum _{j=1}^k|{\mathbf{u}}_j⟩⟨\mathbf{j}|,$$

$${\mathbf{A}}^{\left(i\right)}=\sum _{j=1}^k{\beta }_j^{\left(i\right)}|{\mathbf{u}}_j⟩⟨\mathbf{j}|=\sum _{j=1}^k\frac{{\left({\sigma }_j{\beta }_j^{\left(i-1\right)}\right)}^2+{\lambda }_2}{{\sigma }_j^2{\beta }_j^{\left(i-1\right)}}|{\mathbf{u}}_j⟩⟨\mathbf{j}|$$

$$|{\psi }_0⟩={|0⟩}^a{\left(|0⟩|0⟩\cdots |0⟩\right)}^C{\left(|0⟩|0⟩\cdots |0⟩\right)}^B{\left(|{\psi }_{{\text{A}}^{(i-1)}}⟩\right)}^A.$$

$$|{\psi }_1⟩=\frac{1}{\sqrt{N_1}}{|0⟩}^a\sum _{j=1}^k{\beta }_j^{\left(i-1\right)}{|{\sigma }_j^2⟩}^C{|{\left({\beta }_j^{\left(i-1\right)}\right)}^2⟩}^B|{\mathbf{u}}_j⟩{|{\mathbf{v}}_j⟩}^A.$$

$${\mathbf{U}}_{PE}\left(\mathbf{X}\right)=\left({\mathbf{F}}_{\mathbf{T}}^{†}\otimes \mathbf{I}\right)\left({\sum }_{\tau =0}^{T-1}|\tau ⟩⟨\tau |\otimes e^{i\mathbf{X}\tau t_0/T}\right)\left({\mathbf{H}}^{\otimes t}\otimes \mathbf{I}\right),$$

$$R_f:$$

$$\to {\left[\rho \left(1+\frac{{\lambda }_2}{{\left({\sigma }_j{\beta }_j^{\left(i-1\right)}\right)}^2}\right)|1⟩+\sqrt{1-{\rho }^2{\left(1+\frac{{\lambda }_2}{{\left({\sigma }_j{\beta }_j^{\left(i-1\right)}\right)}^2}\right)}^2}|0⟩\right]}^a{|{\sigma }_j^2⟩}^C{|{\left({\beta }_j^{\left(i-1\right)}\right)}^2⟩}^B,$$

$$|{\psi }_2⟩=$$

$$\otimes \sum _{j=1}^k{\beta }_j^{\left(i-1\right)}{|{\sigma }_j^2⟩}^C{|{\left({\beta }_j^{\left(i-1\right)}\right)}^2⟩}^B|{\mathbf{u}}_j⟩{|{\mathbf{v}}_j⟩}^A.$$

$$|{\psi }_{\mathbf{A}}^{\left(i\right)}⟩=\frac{1}{\sqrt{N_2}}\sum _{j=1}^k\frac{{\left({\sigma }_j{\beta }_j^{\left(i-1\right)}\right)}^2+{\lambda }_2}{{\sigma }_j^2{\beta }_j^{\left(i-1\right)}}|{\mathbf{u}}_j⟩{|{\mathbf{v}}_j⟩}^A.$$

$$t{r}_2\left(|{\psi }_{\mathbf{A}}^{\left(i\right)}⟩⟨{\psi }_{\mathbf{A}}^{\left(i\right)}|\right)=\frac{1}{{\sum }_{k=1}^r{\left({\beta }_j^{\left(i\right)}\right)}^2}{\sum }_{k=1}^r{\left({\beta }_j^{\left(i\right)}\right)}^2|{\mathbf{u}}_k⟩⟨{\mathbf{u}}_k|=\frac{{\text{A}}^{(i)}{{\text{A}}^{(i)}}^{†}}{tr\left({\text{A}}^{(i)}{{\text{A}}^{(i)}}^{†}\right)}.$$

$$y_j=1+\frac{{\lambda }_2}{{\left({\sigma }_j{\beta }_j^{\left(i-1\right)}\right)}^2}.$$

$${\cos }^2\left(2\pi {\omega }_{a_1\cdots a_{i-1}}\right)={\left(\frac{{\alpha }_{a_1{a}_2\cdots a_{i-1}0}}{{\alpha }_{a_1{a}_2\cdots a_{i-1}}}\right)}^2+O\left(poly\left(ϵ\right)\right).$$

$$U_i:$$

$$\stackrel{QRAM}{\to }$$

$$c-S_{{\omega }^{\left(i\right)}}:|{\omega }^{\left(i\right)}⟩|0⟩\to |{\omega }^{\left(i\right)}⟩e^{2\pi i{\omega }^{\left(i\right)}}|0⟩,|{\omega }^{\left(i\right)}⟩|1⟩\to |{\omega }^{\left(i\right)}⟩e^{-2\pi i{\omega }^{\left(i\right)}}|1⟩.$$

$$|{\omega }_l^{\left(i\right)}⟩|0⟩$$

$$|{\omega }_l^{\left(i\right)}⟩|1⟩$$

$$|{\omega }_1^{\left(i\right)}⟩|{\omega }_2^{\left(i\right)}⟩\cdots |{\omega }_p^{\left(i\right)}⟩|0⟩$$

$$|{\omega }_1^{\left(i\right)}⟩|{\omega }_2^{\left(i\right)}⟩\cdots |{\omega }_p^{\left(i\right)}⟩|1⟩$$

$$|{\omega }^{\left(i\right)}⟩|0⟩$$

$$\stackrel{I\otimes H}{\to }$$

$$\stackrel{c-S_{\omega }}{\to }$$

$$\stackrel{I\otimes H}{\to }$$

$$\stackrel{I\otimes V}{\to }$$

$$=$$

$$\theta =arcsin\left(y\right)\approx y+\frac{1}{6}{y}^3+\frac{3}{40}{y}^5+\frac{5}{112}{y}^7+\cdots .$$

$${\tilde{\mathbf{X}}}^T\mathbf{A}$$

$$=\sum _{j=1}^k\sum _{j^{\prime }=1}^k{\sigma }_j{\beta }_{j^{\prime }}^{\left(0\right)}|{\mathbf{v}}_j⟩⟨{\mathbf{u}}_j|{\mathbf{u}}_{j^{\prime }}⟩⟨{\mathbf{j}}^{\prime }|=\sum _{j=1}^k{\sigma }_j{\beta }_j^{\left(0\right)}|{\mathbf{v}}_j⟩⟨\mathbf{j}|$$

$$\mathbf{B}$$

Full text

Quantum machine learning, dimensionality reduction, A-Optimal Projection

I Introduction

Learning low dimensional image representations has gained significant importance in many image processing tasks such as recognition and retrieval Chen et al. (2017); Wang et al. (2017); Song et al. (2017). A range of applications of this problem can be seen in the field of medical imaging such as liver cirrhosis, lung cancer classification and breast cancer diagnosis Song et al. (2017). Another typical example is face recognition which is typically used in security systems or as a commercial identification and marketing tool Song et al. (2018). Recent studies have shown that images are possibly sampled from a low dimensional manifold, however, the visual features, such as color, texture and shape, which are usually extracted for the image representation, are usually of very high dimensionality He et al. (2016). Therefore, a range of techniques have been developed for dimensionality reduction. For instance, principal component analysis (PCA) is guaranteed in terms of the linearly embedded manifold Bishop (2006). Moreover, Isomap, Locally Linear Embedding, and Laplacian Eigenmap are proposed for nonlinear embedded manifold Tenenbaum et al. (2000); Roweis and Saul (2000); Belkin and Niyogi (2001).

Different from all the aforementioned techniques which are not directly related to the regression task, X. He proposed a novel dimensionality reduction algorithm named A-Optimal Projection (AOP) which performs better regression performance in the reduced space He et al. (2016). This approach can be performed under either unsupervised or supervised mode, which is a more widely used algorithm compared to the unsupervised algorithm PCA or the supervised algorithm linear discriminant analysis (LDA) Bishop (2006). Moreover, different from most dimensionality reduction algorithms which are applied as pre-processing of the data, AOP can directly improve the performance of a regression model in the reduced space, therefore, the learned regression model can be as stable as possible.

Time complexity is a significant drawback in classical machine learning algorithms. A range of quantum algorithms have achieved exponentially speed up in machine learning compared with the classical ones Adcock et al. (2015); Wittek (2014). In particular, quantum algorithms for solving the problem of pattern classification and image classification problems were proposed, covering an important area of machine learning Schuld et al. (2016); Duan et al. (2017). Recently, a quantum generative algorithm which is more capable of representing probability distributions was proposed, generating an intriguing link among quantum many-body physics, quantum computational complexity theory and the machine learning frontier Gao et al. (2018). The relationship between feature maps, kernel methods in machine learning and quantum computing was also investigated, and the idea of embedding data into a quantum Hilbert space opens up a promising avenue to quantum machine learning Schuld and Killoran (2019). Moreover, small quantum computers, larger special purpose quantum simulators, annealers, etc., exhibit promising applications in machine learning, and the perspectives on the work of these hardware have also been discussed Biamonte et al. (2017). Quantum machine learning has also been combined with the information security. It was designed to protect private data during performing quantum machine learning, which has potential applications in the big data era Sheng and Zhou (2017).

In the application field of quantum dimensionality reduction, the quantum algorithm for PCA has been proposed for unsupervised mode Lloyd et al. (2014), and the quantum algorithm for LDA has been proposed for supervised mode and classification Cong and Duan (2016). In this paper, we focus on the new dimensionality reduction algorithm AOP which can be used both on unsupervised and supervised model, and propose a quantum algorithm for AOP, which achieves exponentially speedup compared with the classical polylogarithmic in both $n$, the original feature space dimension, and $k$, the reduced feature space dimension.

Our work has two major contributions. First, we present a quantum algorithm for solving the learning process of the AOP algorithm. A reformulation of the original classical AOP algorithm is introduced here which helps the QAOP algorithm be implemented more efficiently. The quantum algorithm is made of iterations, where each iteration mainly consists of phase estimation and a controlled rotation. The reformulated AOP and the partial trace technology can help omit the quantum-classical interactions during the quantum algorithm. Second, we design a detailed quantum circuit for the proposed QAOP algorithm which makes it possible to execute the QAOP algorithm on a universal quantum computer. The circuit for preparing the initial state is presented and the detailed circuit for the controlled rotation is designed. The space and time analysis of the quantum circuit also shows an exponential speedup in the size of the feature space than the classical counterparts.

This paper is arranged as follows: We give a brief overview of the classical AOP algorithm in Sec. II. In Sec. III, the quantum algorithm for AOP algorithm which is used in dimension reduction is presented. In Sec. IV, the overview and detailed quantum circuits for solving QAOP algorithm are designed. Finally we show the conclusions in Sec. V.

II Review of classical A-Optimal Projection

In this section, we briefly review the AOP model and learning algorithm.

The classical AOP dimensionality reduction aims to improve the regression performance in the reduced space which preserves similarities between the data pairs. The AOP dimensionality reduction algorithm returns the directions of projections, and with this result, the data can be projected onto a lower-dimensional subspace which can be directly used for regression problem.

Let ${\bf{X}}=\left({{{\bf{x}}_{1}},\cdots,{{\bf{x}}_{m}}}\right)$ be a $n\times m$ data matrix, where $m$ is the number of data points and $n$ is the number of features. In the graph based dimensionality reduction, we are given a nearest neighbor graph $G$ which represents the geometrical structure of the data manifold. Each vertex of the graph represents a data point. Let ${\bf{S}}\in\mathbb{R}{{}^{m\times m}}$ be the weight matrix of the graph and $N_{k}\left({\bf{x}}\right)$ denote the $k$ nearest neighbors of ${\bf{x}}$. Then a simple example of ${\bf{S}}$ can be defined as follows:

[TABLE]

AOP aims to find a projection matrix ${\bf{A}}\in\mathbb{R}{{}^{n\times k}}$ that maps the the points ${{\bf{x}}_{i}}$ to ${{\bf{y}}_{i}}\in\mathbb{R}{{}^{k}}$ $(i=1,...,m$, and $k\ll n)$, where ${{\bf{y}}_{i}}={{\bf{A}}^{T}}{{\bf{x}}_{i}}$. And using ${{\bf{y}}_{i}}$ to train a linear regression model:

[TABLE]

where $z$ is the observation, ${\bf{\beta}}$ is the weight vector and $\epsilon_{0}$ is an unknown error with Gaussian distribution.

Formally, the objective function of AOP is:

[TABLE]

where ${\lambda_{1}}$ and ${\lambda_{2}}$ are the regularization coefficients which are very small, and ${\bf{L}}=diag\left({{\bf{S1}}}\right)-{\bf{S}}$ is the graph Laplacian (${\bf{1}}$ is a vector of all ones).

To solve the objective function, Ref. He et al. (2016) introduces a variables B and the optimization problem (3) is equivalent to the following:

[TABLE]

where $\widetilde{\bf{X}}={\bf{X}}\Sigma$, and $\Sigma$ is from the cholesky decomposition: ${\bf{I}}+{\lambda_{1}}{\bf{L}}=\Sigma{\Sigma^{T}}$.

It tells us that the optimal ${\bf{A}}$ can be obtained by iteratively computing ${\bf{A}}$ and ${\bf{B}}$. Then the overall procedure of the AOP learning algorithm is depicted as follows:

1. Initialize the matrix ${\bf{A}}$ by computing the PCA of the data ${\bf{X}}$.

2. Compute the matrix ${\bf{B}}$ according to the Eq. (5):

[TABLE]

3. Update the matrix ${\bf{A}}$ according to the Eq. (6), and normalize ${\bf{A}}$ such that ${\left\|{\bf{A}}\right\|_{F}}\leq{\rho_{0}}$.

[TABLE]

where $\rho_{0}$ is used as a constraint parameter to control the size of ${\bf{A}}$.

4. Repeat steps 2 and 3 until convergence.

III Quantum A-Optimal Projection

In this section, we propose the quantum AOP algorithm for dimensionality reduction. We firstly reformulated the original classical AOP algorithm. And with the help of the reformulation, the proposed quantum AOP algorithm can then be implemented more efficiently.

III.1 Reformulation of the AOP algorithm

We reformulated the algorithm in Sec. II in terms of quantum mechanics. Firstly, we adjust the initialization of ${\bf{A}}$ to make it closer to the optimal solution than the original algorithm. Secondly, we combine the steps 2 and 3 into one step and remove the variable B. The advantage of eliminating B is to help avoid quantum-classical transformation during the iteration of the algorithm. Specifically, in one iteration of the original algorithm, a quantum state is needed to be computed and sampled as to reconstruct the matrix B, and then used to update A. In our methods, the elimination of the matrix B can help update the quantum state representing A without sampling and reconstruction in one of the iterations. Finally, by introducing the partial trace technology, quantum-classical interaction can also be omitted between the iterations.

(1) Initialization of ${\bf{A}}$. The original AOP algorithm compute PCA of the data ${\bf{X}}$ to initialize the matrix ${\bf{A}}$. In contrast, we compute PCA of ${\widetilde{\bf{X}}}$ to obtain the initialization ${{\bf{A}}^{\left(0\right)}}$. As shown in Eq. (3), when ${\lambda_{1}}$ and ${\lambda_{2}}$ are set to zero, the objection function (3) is equivalent to the objection function of PCA referring to the data ${\bf{X}}$. And when ${\lambda_{2}}$ is zero, the objection function (3) is equivalent to the objection function of PCA referring to the data ${\widetilde{\bf{X}}}$. It is obvious that the later one is closer to the optimal solution than the former one.

(2) Reformulation of the AOP algorithm. Now turning to the steps 2 and 3 of the AOP algorithm in one of the iterations. Suppose the singular value decomposition of the matrix ${\widetilde{\bf{X}}}$ is ${\widetilde{\bf{X}}}=\sum\nolimits_{j=1}^{r}{{\sigma_{j}}\left|{{{\bf{u}}_{j}}}\right\rangle\left\langle{{{\bf{v}}_{j}}}\right|}$, where $r\leq\min\left({m,n}\right)$ is the rank of ${\widetilde{\bf{X}}}$, and ${\sigma_{k}}\left({{\sigma_{1}}>\cdots>{\sigma_{r}}>0}\right)$ are the singular values of ${\widetilde{\bf{X}}}$, with ${{\bf{u}}_{j}}$ and ${{\bf{v}}_{j}}$ being the left and right singular vectors. Obviously, we have ${\widetilde{\bf{X}}^{T}}=\sum\nolimits_{j=1}^{r}{{\sigma_{j}}\left|{{{\bf{v}}_{j}}}\right\rangle\left\langle{{{\bf{u}}_{j}}}\right|}$ and $\widetilde{\bf{X}}{\widetilde{\bf{X}}^{T}}=\sum\nolimits_{j=1}^{r}{\sigma_{j}^{2}\left|{{{\bf{u}}_{j}}}\right\rangle\left\langle{{{\bf{u}}_{j}}}\right|}$.

As ${{\bf{A}}^{\left(0\right)}}$ is the PCA of ${\widetilde{\bf{X}}}$, we have

[TABLE]

where $k$ is the rank of ${\bf{A}}$, and $\left|{\bf{j}}\right\rangle$’s are the basis states. Now we have the theorem 1 (and the proof is shown in Appendix A.).

Theorem 1:

• Given the matrix ${\widetilde{\bf{X}}}=\sum\nolimits_{j=1}^{r}{{\sigma_{j}}\left|{{{\bf{u}}_{j}}}\right\rangle\left\langle{{{\bf{v}}_{j}}}\right|}$, the $i$-th iteration of the AOP algorithm outputs the matrix ${{\bf{A}}^{\left(i\right)}}$:*

[TABLE]

*where $1\leq i\leq s$ and ${\beta_{j}^{\left(0\right)}}=1$ for all $j$’s. *

According to the theorem 1, the reformulated AOP algorithm is presented as follows:

1. Initialize the matrix ${{\bf{A}}^{\left(0\right)}}$ by computing the PCA of ${\widetilde{\bf{X}}}$ according to Eq. (7).

2. Update the matrix ${\bf{A}}$ according to Eq. (8).

3. Repeat step 2 until convergence.

The modeling of the AOP algorithm is shown in Fig. 1.

III.2 QAOP algorithm

The overall procedure of our QAOP algorithm is then proposed as follows.

Algorithm. $\textbf{A}^{(s)}=QAOP\left({\widetilde{\bf{X}}},{\lambda_{2}}\right)$.

1. Initialize $i=1$ , apply quantum PCA algorithm of ${\widetilde{\bf{X}}}$ to compute $\textbf{A}^{(0)}$ Lloyd et al. (2014).

2. Perform one iteration $Q$ of the QAOP algorithm to compute $\textbf{A}^{(i)}$, i.e. $\textbf{A}^{(i)}=Q\left({\widetilde{\bf{X}}},\textbf{A}^{(i-1)}\right)$.

3. Set $i=i+1$, and repeat step 2 until the number of iterations $i=s$; now the projection matrix $\textbf{A}^{(s)}$ can be obtained.

The quantum algorithm for $i$-th iteration $Q$ of the QAOP algorithm in Step 2 is presented as follows.

(1) Prepare four quantum registers in the state

[TABLE]

where the superscript $a$ represents the ancilla qubit, the superscripts ${C},{B},A$ represent the register ${C},{B}$ and $A$, respectively.

(2) Perform the unitary operation ${{\bf{U}}_{PE}}\left({\widetilde{\bf{X}}}{\widetilde{\bf{X}}}^{\dagger}\right)$ and ${{\bf{U}}_{PE}}\left({\textbf{A}^{(i-1)}}{\textbf{A}^{(i-1)}}^{\dagger}\right)$ on the state, then we have the state

[TABLE]

Here ${{\bf{U}}_{PE}}$ represents the unitary matrix for phase estimation which we fully characterized in Duan et al. (2017):

[TABLE]

where ${\bf{F}}_{{\bf{T}}}^{\dagger}$ is the inverse quantum Fourier transform and ${\sum\nolimits_{\tau=0}^{T-1}{\left|\tau\right\rangle{{\left\langle\tau\right|}^{C}}\otimes{e^{i{\bf{A}}\tau{{{t_{0}}}\mathord{\left/{\vphantom{{{t_{0}}}T}}\right.\kern-1.2pt}T}}}}}$ is the conditional Hamiltonian evolution Harrow et al. (2009).

(3) Apply a controlled rotation ${\bf{R}}_{f}$ to the ancilla qubit, controlled by both the register ${C}$ and ${B}$. ${\bf{R}}_{f}$ is defined as follows:

[TABLE]

where $\rho<1$ is a parameter used for normalization of the quantum state.

This rotation transforms the state to

[TABLE]

(4) Uncompute the registers $C$, $B$ and $A$, remove the register $C$ and $B$, and measure the ancilla qubit to be $\left|1\right\rangle$. Then, we have the state proportional to

[TABLE]

Here, we can construct the matrix ${\textbf{A}^{(i)}}{\textbf{A}^{(i)}}^{\dagger}$ for the phase estimation in the next iteration by taking a partial trace of $\left|{\psi_{\bf{A}}^{\left(i\right)}}\right\rangle\left\langle{\psi_{\bf{A}}^{\left(i\right)}}\right|$. Note that the eigenvectors of ${\textbf{A}^{(i)}}{\textbf{A}^{(i)}}^{\dagger}$ are ${{\bf{u}}_{j}}$ and the corresponding eigenvalues are ${\left({\beta_{j}^{\left({i}\right)}}\right)}^{2}$. Then the density matrix that represents ${\textbf{A}^{(i)}}{\textbf{A}^{(i)}}^{\dagger}$ can be obtained Schuld et al. (2016) :

[TABLE]

IV QAOP circuit

In this section, we study the QAOP algorithm in terms of the quantum circuit model. The quantum circuits provide the possibility to implement the quantum algorithm on a universal quantum computer. First, we present the overview model of the QAOP circuit. Second, we study in depth the realization of the initial state preparation and the controlled rotation in terms of the quantum circuit model. Finally, the space and time resources required for the quantum circuit are analyzed.

The overview of the circuit for solving QAOP is shown in Fig. 2. It provides one model to implement the QAOP algorithm. Take the $i$-th iteration of the QAOP algorithm for example. It can be divided into three major steps: (1) Phase estimation: as the eigenspace of the unitary ${e^{-i\widetilde{\bf{X}}{{\widetilde{\bf{X}}}^{T}}{t_{0}}}}$ and ${e^{-i{{\bf{A}}^{\left({i-1}\right)}}{{\left({{{\bf{A}}^{\left({i-1}\right)}}}\right)}^{T}}{t_{0}}}}$ are both spanned by the eigenvectors ${\left|{{{\bf{u}}_{j}}}\right\rangle}$, they can be both implemented on the input state ${\left|{\psi_{\bf{A}}^{\left({i-1}\right)}}\right\rangle}$. And the Hadamard gates and the inverse QFT of ${{\bf{U}}_{PE}}\left({\widetilde{\bf{X}}}{\widetilde{\bf{X}}}^{\dagger}\right)$ and ${{\bf{U}}_{PE}}\left({\textbf{A}^{(i-1)}}{\textbf{A}^{(i-1)}}^{\dagger}\right)$ can be implemented in parallel. (2) Controlled rotation: it consists of ${U_{\beta,\sigma}}$ and $c-R_{y}$. Firstly the operation ${U_{\beta,\sigma}}$ computes the function ${y_{j}}$ of the output eigenvalues of ${\widetilde{\bf{X}}}{\widetilde{\bf{X}}}^{\dagger}$ and ${\textbf{A}^{(i-1)}}{\textbf{A}^{(i-1)}}^{\dagger}$ as follows:

[TABLE]

And then the operation $c-R_{y}$ extracts the value of $y_{j}$ in the basic states of the register $L$ to the amplitude of the ancilla qubit. (3) Uncomputing: undo the Reg. $B,C$ and $L$, and measure the top ancilla qubit. If the result returns to 1, then the Reg. $A$ of the quantum system collapses to the output state ${\left|{\psi_{\bf{A}}^{\left({i}\right)}}\right\rangle}$, which is also the input state of the $\left(i+1\right)$-th iteration.

We now deal with the detailed QAOP circuit. In the following, we mainly investigate the quantum circuits for the initial state preparation and the controlled rotation in one iteration of the QAOP circuit.

IV.1 State preparation

At the very beginning, we present a detailed quantum circuit for preparing the initial state of the QAOP algorithm. Suppose each element of ${\textbf{A}^{(0)}}\in\mathbb{R}{{}^{n\times k}}$ is given, and its corresponding quantum state ${\left|{{\psi_{\textbf{A}^{(0)}}}}\right\rangle}=\left|{{a_{1}}{a_{2}}\cdots{a_{q}}}\right\rangle$ is a $q$-qubit quantum state, where $q=O\left[\log_{2}\left(nk\right)\right]$. Following the approach in Mosca and Kaye (2001), the quantum circuit for the initial state preparation can be shown in Fig. 3. Here, we introduce quantum random access memory (QRAM) to omit the register $\left|{\bar{\psi}}\right\rangle$ in Mosca and Kaye (2001).

In Fig. 3, an register of $p=O\left(\log_{2}\epsilon^{-1}\right)$ qubits is used for storing the ${\omega^{\left(i\right)}}$, where ${\omega^{\left(i\right)}}$ is short for $\omega_{{a_{1}}\cdots{a_{i-1}}}\left(i=1,2,3,...,q\right)$ satisfying:

[TABLE]

All the $\omega_{{a_{1}}\cdots{a_{i-1}}}$ can be computed classically and they are supposed to be stored in the QRAM. Given the index ${{a_{1}}{a_{2}}\cdots{a_{i-1}}}$, define $j=h\left({{a_{1}}{a_{2}}\cdots{a_{i-1}}}\right)$ being the address where the data ${\omega_{{a_{1}}{a_{2}}\cdots{a_{i-1}}}}$ stores, where $h\left(\cdot\right)$ is a hash function mapping ${{a_{1}}{a_{2}}\cdots{a_{i-1}}}$ to $j$. Then define the unitary operation $U_{i}$ which implements the QRAM readout operation Giovannetti et al. (2008):

[TABLE]

specifically, $U_{i}$ outputs the content ${\omega_{{a_{1}}{a_{2}}\cdots{a_{i-1}}}}$ of the $j$-th memory cell in QRAM. Ref. Hong et al. (2012) shows that this procedure can be implemented in time $O\left(p\right)$.

And the number of ${U_{i}}$ in the circuit for state preparation is $q$, so the total memory calls of QRAM is $O\left(pq\right)$. Therefore, an inverse in $pq$ error rate suffices to achieve an overall constant error per QRAM look-up Arunachalam et al. (2015).

Moreover, we make further study on $c-{S_{\omega^{\left(i\right)}}}$ which is defined in Ref. Mosca and Kaye (2001):

[TABLE]

Define {R_{l}}=\left({\begin{array}[]{*{20}{c}}{{e^{{{2\pi i}\mathord{\left/{\vphantom{{2\pi i}{{2^{l}}}}}\right.\kern-1.2pt}{{2^{l}}}}}}}&0\\ 0&{{e^{{{-2\pi i}\mathord{\left/{\vphantom{{-2\pi i}{{2^{l}}}}}\right.\kern-1.2pt}{{2^{l}}}}}}}\end{array}}\right), then the controlled unitary $c-R_{l}$ implements the following transformation:

[TABLE]

where $\omega_{l}^{\left(i\right)}$ is the $l$-th binary bit of $\omega^{\left(i\right)}$, specifically, ${\omega^{\left(i\right)}}={2^{-1}}\omega_{1}^{\left(i\right)}+{2^{-2}}\omega_{2}^{\left(i\right)}+\cdots+{2^{-p}}\omega_{p}^{\left(i\right)}=0.\omega_{1}^{\left(i\right)}\omega_{2}^{\left(i\right)}\cdots\omega_{p}^{\left(i\right)}$. Now $\prod\limits_{l=1}^{p}{\left({c-{R_{l}}}\right)}$ achieves the function of $c-{S_{\omega^{\left(i\right)}}}$:

[TABLE]

Therefore, the quantum circuit for $c-{S_{\omega^{\left(i\right)}}}$ can be implemented as shown in Fig. 4.

Now we can simply infer that the number of qubits needed for preparing the initial quantum state ${\left|{{\psi_{\textbf{A}^{(0)}}}}\right\rangle}$ is $O\left(p+q\right)$, and the number of gates required is $O\left(pq\right)$.

In summary, with V=\left({\begin{array}[]{*{20}{c}}1&0\\ 0&{-\iota}\end{array}}\right), the unitary $\left({I\otimes V}\right)\left({I\otimes H}\right)\left({c-{S_{\omega}}}\right)\left({I\otimes H}\right)$ implements the transformation:

[TABLE]

IV.2 The controlled rotation

The controlled rotation mainly involves the computation of $U_{\beta,\sigma}$ and $c-R_{y}$. In the stage of $U_{\beta,\sigma}$, Newton’s method is introduced for computing $y_{j}=y\left({\sigma_{j}^{2},\beta_{j}^{2}}\right)=\frac{\rho\left({\sigma_{j}^{2}\beta_{j}^{2}+{\lambda_{2}}}\right)}{{\sigma_{j}^{2}\beta_{j}^{2}}}=\rho+\frac{{\rho{\lambda_{2}}}}{{\sigma_{j}^{2}\beta_{j}^{2}}}$, where $\rho<1$ is used for normalization of the quantum state. The value of $y_{j}$ are computed out and stored in the basis state of the register $L$, and the number of qubits for storing $y_{j}$ is $d=O\left({{{\log}_{2}}\kappa}\right)$. In the stage of $c-R_{y}$, $y_{j}$ is used as controlled qubit controlling the top ancilla qubit in Fig. 2.

(1) For $U_{\beta,\sigma}$, let ${z_{j}}=z\left({\sigma_{j}^{2}\beta_{j}^{2}}\right)={1\mathord{\left/{\vphantom{1{\left({\sigma_{j}^{2}\beta_{j}^{2}}\right)}}}\right.\kern-1.2pt}{\left({\sigma_{j}^{2}\beta_{j}^{2}}\right)}}$, then we have ${y_{j}}=\rho+\rho{\lambda_{2}}z_{j}^{\left({s^{\prime}}\right)}$. Here we use Newton iteration to approximate ${1\mathord{\left/{\vphantom{1{\left({\sigma_{j}^{2}\beta_{j}^{2}}\right)}}}\right.\kern-1.2pt}{\left({\sigma_{j}^{2}\beta_{j}^{2}}\right)}}$, for ${\sigma_{j}^{2}\beta_{j}^{2}}>1$. The quantum circuit for computing the initial approximation ${z_{j}^{\left(0\right)}}$ can be seen in Cao et al. (2012). Applying the Newton method to $f\left({{z_{j}}}\right)={1\mathord{\left/{\vphantom{1{{z_{j}}}}}\right.\kern-1.2pt}{{z_{j}}}}-\sigma_{j}^{2}\beta_{j}^{2}$, we can get the Newton iteration function: $z_{j}^{\left({i+1}\right)}=g\left({z_{j}^{\left(i\right)}}\right)=z_{j}^{\left(i\right)}-\frac{{f\left({z_{j}^{\left(i\right)}}\right)}}{{f^{\prime}\left({z_{j}^{\left(i\right)}}\right)}}=-\sigma_{j}^{2}\beta_{j}^{2}{\left({z_{j}^{\left(i\right)}}\right)^{2}}+2z_{j}^{\left(i\right)}$. Then the detailed circuit for one iteration of Newton’s method is presented in Fig. 5. The number of qubits needed in the ancilla registers is $O(d)$. The number of fundamental quantum operations for implementing addition and multiplication is $O\left[poly\left(d\right)\right]$, where the degree of the polynomial is no more than $3$ Vedral et al. (1996); Parent et al. (2018), and the number of gates for implementing shift is $O(d)$.

The quantum circuit for ${y_{j}}=\rho+\rho{\lambda_{2}}z_{j}^{\left({s^{\prime}}\right)}$ can be simply realized as shown in Fig. 6. The circuit can be simply realized with the quantum circuits for addition and multiplication. Therefore, the number of qubits and gates needed in the circuit are $O(d)$ and $O\left[poly\left(d\right)\right]$, respectively.

To sum up, the overall quantum circuit for $U_{\beta,\sigma}$ can be designed as shown in Fig. 7. And we can simply infer that the number of qubits needed in these circuits is $O(d+b)$. Let the number of Newton iteration be $s^{\prime}$, then the number of gates required in Fig. 7 is $O\left[s^{\prime}poly\left(d\right)\right]$.

Now we analyze the error caused by Newton’s iteration. Similar to the error analysis in Duan et al. (2018); Cao et al. (2012), the error consists of two parts. One is error $e_{s^{\prime}}$ caused by the Newton’s iteration, the other is the roundoff error ${{\hat{e}}_{s^{\prime}}}$ caused by truncating the result of one iteration to $d$ qubits of accuracy.

According to the Newton iteration function, we have $g\left({z_{j}^{\left(i\right)}}\right)-\frac{1}{{\sigma_{j}^{2}\beta_{j}^{2}}}=-\sigma_{j}^{2}\beta_{j}^{2}{\left({z_{j}^{\left(i\right)}-\frac{1}{{\sigma_{j}^{2}\beta_{j}^{2}}}}\right)^{2}}$. Then the error $e_{s^{\prime}}$ satisfies ${e_{s^{\prime}}}:=\left|{z_{j}^{\left({s^{\prime}}\right)}-\frac{1}{{\sigma_{j}^{2}\beta_{j}^{2}}}}\right|=\sigma_{j}^{2}\beta_{j}^{2}e_{{s^{\prime}}-1}^{2}=\frac{1}{{\sigma_{j}^{2}\beta_{j}^{2}}}{\left({\sigma_{j}^{2}\beta_{j}^{2}{e_{0}}}\right)^{{2^{s^{\prime}}}}}$. Following the approach in Cao et al. (2012), the initial error ${e_{0}}$ satisfies ${\sigma_{j}^{2}\beta_{j}^{2}}{e_{0}}<1/2$, then for error ${\varepsilon_{N}}$ we have ${2^{-{2^{s^{\prime}}}}}\leq{\varepsilon_{N}}$, which implies ${s^{\prime}}\geq\left\lceil{{{\log}_{2}}{{\log}_{2}}\varepsilon_{N}^{-1}}\right\rceil$, where ${\varepsilon_{N}}$ denotes the desired error of Newton iteration without considering the truncation error. We can also follow the result in Cao et al. (2012) that the truncation error ${{\hat{e}}_{s^{\prime}}}$ satisfies ${{\hat{e}}_{s^{\prime}}}:=\left|{{{\hat{z}}^{\left({s^{\prime}}\right)}}-{z^{\left({s^{\prime}}\right)}}}\right|\leq{s^{\prime}}2^{-d}$.

In short, with the number of the iteration steps being ${s^{\prime}}=O\left({{{\log}_{2}}d}\right)$, the error caused by the unitary ${\bf{U}}_{\sigma,\tau}$ is $\left|{z_{j}^{\left({s^{\prime}}\right)}-\frac{1}{{\sigma_{j}^{2}\beta_{j}^{2}}}}\right|\leq{\varepsilon_{N}}+{s^{\prime}}2^{-d}$, where ${\varepsilon_{N}}\geq{2^{-{2^{s^{\prime}}}}}$.

(2) For $c-R_{y}$, in order to make the output quantum state accurate, the rotation angle $\theta$ of $R_{y}$ satisfies $\theta=arc\sin\left(y\right)$. Since $arcsin$ has a convergent Taylor series, we can approximate

[TABLE]

Then the quantum circuit for computing the rotation angle of ${c-R_{y}}$ can be implemented as shown in Fig. 8. This circuit also only consists of the operations for addition and multiplication, therefore the number of qubits and gates required are $O(d)$ and $O\left[poly\left(d\right)\right]$ respectively.

The quantum circuit for $c-R_{y}$ is shown in Fig. 9, where $\theta_{1},\cdots,\theta_{d}$ are the binary bits of the output $\theta$ in Fig. 8. Obviously, the space and time complexity are both $O(d)$.

IV.3 Complexity

We firstly analyze the space and time resources used in phase estimation. Let the efficient condition number of ${\widetilde{\bf{X}}}{\widetilde{\bf{X}}}^{\dagger}$ be $\kappa$. As ${\textbf{A}^{(0)}}$ is induced by PCA of ${\widetilde{\bf{X}}}$, the condition number of ${\textbf{A}^{(0)}}{\textbf{A}^{(0)}}^{\dagger}$ is not greater than $\kappa$, and so as the ${\textbf{A}^{(i)}}{\textbf{A}^{(i)}}^{\dagger}$ for $i=1,\cdots,s-1$. Therefore, ${\sigma_{j}}^{2}\in\left[{1/\kappa,1}\right]$ and ${\beta_{j}}^{2}\in\left[{1/\kappa,1}\right]$, and the numbers of qubits to store $\sigma_{j}^{2}$ in Reg. $B$ and $\beta_{j}^{2}$ in Reg. $C$ are both $b=O\left({{{\log}_{2}}\kappa}\right)$. Moreover, we can learn from Nielsen and Chuang (2010) that $O\left({{b^{2}}}\right)$ operations and two calls to the controlled-unitary black boxes are needed in the stage of phase estimation.

We now analyze the space and time complexity of the whole QAOP circuit. In one iteration $Q$ of the QAOP circuit, the number of qubits required for preparing the initial quantum state ${\left|{{\psi_{\textbf{A}^{(0)}}}}\right\rangle}$ is $O\left(p+q\right)=O\left[\log_{2}\left(nk/\epsilon\right)\right]$. And phase estimation requires $O\left({{{\log}_{2}}\kappa}\right)$ qubits, where the condition number $\kappa$ is usually taken as $\kappa=O\left({1/\epsilon}\right)$. Therefore, the space complexity of phase estimation is $O\left(\log_{2}\epsilon^{-1}\right)$. Taking the ancilla qubits in the controlled rotation into account, the number of qubits in this stage is $O(d+b)=O\left({{{\log}_{2}}\kappa}\right)=O\left(\log_{2}\epsilon^{-1}\right)$. The number of qubits will not increase with the number of iterations of QAOP. To sum up, the total number of qubits required in the quantum circuit is $O\left[\log_{2}\left(nk/\epsilon\right)\right]$.

Now turning to the time consumption. In one iteration $Q$, the number of gates for the initial state preparing stage is $O\left(pq\right)=O\left[\log_{2}(nk)\log_{2}\epsilon^{-1}\right]$. And phase estimation requires $O\left(b^{2}\right)=O\left[\left({{{\log}_{2}}\kappa}\right)^{2}\right]=O\left[\left({{{\log}_{2}}\epsilon^{-1}}\right)^{2}\right]$ operations and two calls to the controlled-unitary black boxes. The number of quantum gates in the controlled rotation is $O\left[s^{\prime}poly\left(d\right)\right]=O\left[s^{\prime}{{poly}}\left({{\log}_{2}}\kappa\right)\right]=O\left[s^{\prime}{{poly}}\left({{\log}_{2}}\epsilon^{-1}\right)\right]$. The number of iteration $Q$ is $s$, therefore the total time complexity of the QAOP circuit is $O[ss^{\prime}{\log_{2}(nk)}{{poly}}\left({{\log}_{2}}\epsilon^{-1}\right)]$.

In summary, $O\left[{{{\log}_{2}}\left({nk}/{\epsilon}\right)}\right]$ space and $O[{\log_{2}(nk)}{{poly}}\left({{\log}_{2}}\epsilon^{-1}\right)]$ elementary operations allow us to implement the QAOP circuit with fidelity at least $1-\epsilon$, when the number of iterations of the QAOP algorithm $s$ and the number of Newton’s iteration $s^{\prime}$ are both small constant numbers.

V Conclusions

We have shown that the proposed quantum algorithm QAOP can be used to speed up the learning process of an important dimensionality reduction algorithm in pattern recognition and machine learning. We firstly reformulated the original AOP algorithm, therefore the quantum-classical interactions during the quantum algorithm can be omitted. We then proposed the QAOP algorithm and investigated the quantum circuits for solving the QAOP algorithm. The detailed quantum circuits for preparing the input quantum state and the circuits for solving the controlled rotation are presented. The space and time complexity of the quantum circuit show that the number of the qubits and gates required are $O\left[{{{\log}_{2}}\left({nk}/{\epsilon}\right)}\right]$ and $O[{\log_{2}(nk)}{{poly}}\left({{\log}_{2}}\epsilon^{-1}\right)]$, respectively. The result shows that the QAOP algorithm and QAOP circuit for dimensionality reduction may motivate to conduct new investigations in quantum machine learning.

Acknowledgements.

This work was supported by the Funding of National Natural Science Foundation of China (Grants No. 61571226 and No. 61701229), Natural Science Foundation of Jiangsu Province, China (Grant No. BK20170802), China Postdoctoral Science Foundation funded Project (Grants No. 2018M630557 and No. 2018T110499), Jiangsu Planned Projects for Postdoctoral Research Funds (Grant No. 1701139B), and the Nanjing University of Aeronautics and Astronautics PhD short-term visiting scholar project. The authors also acknowledge Michele Mosca for inspiring discussions.

Appendix A Proof of theorem 1

Proof: According to Eq. (7), we have the initialization of ${{\bf{A}}^{\left(0\right)}}=\sum\limits_{j=1}^{k}{{\beta_{j}^{\left(0\right)}}\left|{{{\bf{u}}_{j}}}\right\rangle\left\langle{\bf{j}}\right|}$, where ${\beta_{j}^{\left(0\right)}}=1$ for all $j$’s. Obviously, $\left\langle{{{{\bf{u}}_{j}}}}{\left|{\vphantom{{{{\bf{u}}_{j}}}{{{\bf{u}}_{j^{\prime}}}}}}\right.\kern-1.2pt}{{{{\bf{u}}_{j^{\prime}}}}}\right\rangle=0$ when $j\neq j^{\prime}$ and $\left\langle{{{{\bf{u}}_{j}}}}{\left|{\vphantom{{{{\bf{u}}_{j}}}{{{\bf{u}}_{j^{\prime}}}}}}\right.\kern-1.2pt}{{{{\bf{u}}_{j^{\prime}}}}}\right\rangle=1$ when $j=j^{\prime}$. In the following, we simplify ${{\bf{A}}^{\left(0\right)}}$ as ${\bf{A}}$. Therefore, we have

[TABLE]

and ${\widetilde{\bf{X}}^{T}}{\bf{A}}{{\bf{A}}^{T}}\widetilde{\bf{X}}=\sum\limits_{j=1}^{k}{{\left(\sigma_{j}{\beta_{j}^{\left(0\right)}}\right)^{2}}\left|{{{\bf{v}}_{j}}}\right\rangle\left\langle{{{\bf{v}}_{j}}}\right|}$. According to the Eq. (5), we have:

[TABLE]

Similarly, we have

[TABLE]

and $\widetilde{\bf{X}}{\bf{B}}{{\bf{B}}^{T}}{\widetilde{\bf{X}}^{T}}=\sum\limits_{j=1}^{k}{{\left(\frac{{{\sigma_{j}^{2}{\beta_{j}^{\left(0\right)}}}}}{{{\left(\sigma_{j}{\beta_{j}^{\left(0\right)}}\right)^{2}}+{\lambda_{2}}}}\right)^{2}}\left|{{{\bf{u}}_{j}}}\right\rangle\left\langle{{{\bf{u}}_{j}}}\right|}$.

According to the Eq. (6), we have:

[TABLE]

Eq. (27) shows that after one iteration of the algorithm, only the singular values of ${\bf{A}}$ are updated while the singular vectors stay the same. Therefore, for $i$-th iteration of the algorithm, let ${{\bf{A}}^{\left(i\right)}}=\sum\limits_{j=1}^{k}{{\beta_{j}^{\left(i\right)}}\left|{{{\bf{u}}_{j}}}\right\rangle\left\langle{\bf{j}}\right|}$, where ${\beta_{j}^{\left(0\right)}}=1$ for all $j$’s according to Eq. (7). We can easily have $\beta_{j}^{\left(i\right)}=\frac{{{{\left({{\sigma_{j}}\beta_{j}^{\left({i-1}\right)}}\right)}^{2}}+{\lambda_{2}}}}{{\sigma_{j}^{2}\beta_{j}^{\left({i-1}\right)}}}$.

∎