Quantum Gaussian process regression

The specific expression is $\left|{{{\vec{y}}}}\right\rangle_{\rm{1}}{\rm{=}}\sum\limits_{j=0}^{M{\rm{-}}1}{\beta_{j}\left|{\vec{u}_{j}}\right\rangle_{1}}$ and $\left|{{{\vec{k}}}_{*}}\right\rangle_{\rm{1}}{\rm{=}}\sum\limits_{j=0}^{M{\rm{-}}1}{\alpha_{j}\left|{\vec{u}_{j}}\right\rangle_{1}}$. Ref. [28] pointed that the procedure of preparing quantum state $\left|{{{\vec{y}}}}\right\rangle_{\rm{1}}$ and $\left|{{{\vec{k}}}_{*}}\right\rangle_{\rm{1}}$ is as follows. Firstly, supposing that the quantum oracles $O_{\vec{y}}$ and $O_{\vec{k}_{*}}$ are provided, which can efficiently access the entries of $\vec{y}$ and $\vec{k}_{*}$ in time $O\left({{\text{poly}}\log M}\right)$ from quantum random access memory (QRAM) and act as $O_{\vec{y}}:\left|j\right\rangle_{1}\left|0\right\rangle_{a1}\mapsto\left|j\right\rangle_{1}\left|{y_{j}}\right\rangle_{a1}$ and $O_{\vec{k}_{*}}:\left|j\right\rangle_{1}\left|0\right\rangle_{b1}\mapsto\left|j\right\rangle_{1}\left|{k_{*j}}\right\rangle_{b1}$. Secondly, giving initial state $\left|0\right\rangle_{1}\left|0\right\rangle_{a1}$ and performing the quantum Fourier transform on the first register, which gets $\sum\limits_{j=0}^{M-1}{\frac{1}{{\sqrt{M}}}\left|j\right\rangle_{1}}\left|0\right\rangle_{a1}$. The quantum Fourier transform on an orthonormal basis $\left|0\right\rangle,...,\left|{N-1}\right\rangle$ is defined to be a linear operator with following action on the basis states,

Thirdly, applying the oracle $O_{\vec{y}}$ on $\sum\limits_{j=0}^{M-1}{\frac{1}{{\sqrt{M}}}\left|j\right\rangle_{1}}\left|0\right\rangle_{a1}$, which obtains $\sum\limits_{j=0}^{M-1}{\frac{1}{{\sqrt{M}}}\left|j\right\rangle_{1}}\left|{y_{j}}\right\rangle_{a1}\textbf{}$. Then, appending one qubit and performing conditioned rotation on the register a1, which gains state

Finally, performing inverse oracle operation and measuring the last register in the basis $\left\{{\left|0\right\rangle,\left|1\right\rangle}\right\}$, the remainder particles are in the state $\frac{{\sum\limits_{j=0}^{M-1}{y_{j}\left|j\right\rangle_{1}}}}{{\left\|{\vec{y}}\right\|_{\max}}}$, which has probability $p_{{y}}=\frac{{\sum\limits_{j=0}^{M-1}{y_{j}^{2}}}}{{M\left\|{{{\vec{y}}}}\right\|_{\max}^{2}}}{\rm{=}}\Omega\left(1\right)$ ($\vec{y}$ and ${{\vec{k}}}_{*}$ are balanced). Thus, we need to perform $O\left(1\right)$ to get $\left|{\vec{y}}\right\rangle$. Moreover, $\left\|{\vec{y}}\right\|^{2}=\sum\limits_{j=0}^{M-1}{y_{j}^{2}}=MP_{y}\left\|{{{\vec{y}}}}\right\|_{\max}^{2}$ can also be estimated. The total runtime of above operation is $O\left({{\text{poly}}\log M}\right)$. In a similar way, $\left|{{{\vec{k}}}_{*}}\right\rangle$ and $\left\|{{{\vec{k}}}_{*}}\right\|^{2}$ can be obtained in time $O\left({{\text{poly}}\log M}\right)$.

Step B1 Assuming that covariance matrix $K$ is prepared (The detailed preparation method will be described in section E). The matrix $K$ is Hermitian but not necessarily sparse. Thus, the unitary linear decomposition (LCU) approach [24] is adopted to take the exponentiation of matrix. LCU means under the assumptions that the quantum oracle $O_{\vec{k}}\left|0\right\rangle^{\otimes\left\lceil{\log M}\right\rceil}=\sum\nolimits_{i=0}^{M-1}{\sqrt{k_{i}}}\left|i\right\rangle$, which can be efficiently implemented in time $O\left({{\rm{poly}}\log M}\right)$ and $\sum\nolimits_{i=0}^{M-1}{k_{i}=1}\textbf{}$ is provided. Zhou and Wang [22] proposed that $e^{-iKt}$ can be simulated within spectral-norm error $\varepsilon$ in time $O\left({\frac{{t{\text{poly}}\log M\log\left({\frac{t}{\varepsilon}}\right)}}{{\log\log\left({\frac{t}{\varepsilon}}\right)}}}\right)$. Noted that the algorithm decomposes $K$ into $M$ linear combinations of unitary operators that can be efficiently implemented, that is $K=\sum\nolimits_{j=0}^{M-1}{k_{j}V_{j}}$, where $V_{j}=\sum\nolimits_{l=0}^{M-1}{\left|{\left({l-j+1}\right)\bmod\left|M\right\rangle\left\langle l\right|}\right|},j=1,2,...,M$ can be implemented using $O\left({\log M}\right)$ one- or two-qubit gates.

Step B2 Performing the phase estimation [32] on $\left|{{{\vec{k}}}_{*}}\right\rangle_{1}=\sum\limits_{j=0}^{M{\rm{-}}1}{\alpha_{j}\left|{u_{j}}\right\rangle_{1}}$, we can obtain

Step B3 Adding an ancilla qubit $\left|0\right\rangle_{3}$ and performing controlled rotation on it controlled on the second register, the state becomes

Step B4 Performing inverse phase estimation and discarding the second register, we can get

The above operations are all unitary operations, thus, above processes can be denoted with $U$.

Step C1 Preparing quantum initial state $\left|0\right\rangle_{1}\left|0\right\rangle_{2}\left({\left|0\right\rangle_{3}\left|{{{\vec{k}}}_{*}}\right\rangle_{4}+\left|1\right\rangle_{3}\left|{{{\vec{y}}}}\right\rangle_{4}}\right)\left|0\right\rangle_{5}$. When the third register is $\left|0\right\rangle_{3}$, unitary operation $U$ on the fourth and fifth register is executed. Otherwise, pauli-X on the fifth register is performed. Thus, the state of system is

Step C2 Performing the Hadamard transform on the first register. In the meanwhile, executing pauli-X on the second register. Then, the Hadamard transform on the second register is carried out. In this way, the state becomes

Step C3 Doing nothing when the first register is $\left|0\right\rangle_{1}$. Otherwise, performing swap operation [25] on the second and the third register, which obtains

Step C4 Measuring the first register of Eq. (17) with operator $M=\frac{1}{2}\left({\left|0\right\rangle-\left|1\right\rangle}\right)\left({\left\langle 0\right|-\left\langle 1\right|}\right)$, where the inner product of Eq. (14) and $\left|{\vec{y}}\right\rangle_{4}\left|1\right\rangle_{5}$ is $\sum\limits_{j=0}^{M-1}{\frac{{c\alpha_{j}\beta_{j}}}{{\lambda_{j}+\sigma_{M}^{2}}}}$. Thus, the probability of measurement is $p=\frac{1}{2}\left({1+\sum\limits_{j=0}^{M-1}{\frac{{c\alpha_{j}\beta_{j}}}{{\lambda_{j}+\sigma_{M}^{2}}}}}\right)$. Here, $c$ is a constant. Through measurement, the specific value of $\sum\limits_{j=0}^{M-1}{\frac{{c\alpha_{j}\beta_{j}}}{{\lambda_{j}+\sigma_{M}^{2}}}}$ is attained. Therefore, predicting mean value in quantum context is achieved.

Step D1 Preparing quantum initial state $\left|0\right\rangle_{1}\left|0\right\rangle_{2}\left({\left|0\right\rangle_{3}\left|{{{\vec{k}}}_{*}}\right\rangle_{4}+\left|1\right\rangle_{3}\left|{{{\vec{k}}}_{*}}\right\rangle_{4}}\right)\left|0\right\rangle_{5}$. When the third register is in the state $\left|0\right\rangle_{3}$, unitary operation $U$ on the fourth register is performed. Otherwise, pauli-X operation is performed on the fifth register. Thus, we can get $\left|0\right\rangle_{1}\left|0\right\rangle_{2}\left({\left|0\right\rangle_{3}\left|{\varphi_{3}}\right\rangle_{45}+\left|1\right\rangle_{3}\left|{{{\vec{k}}}_{*}}\right\rangle_{4}\left|1\right\rangle_{5}}\right)$.

Step D2 Performing the Hadamard transform on the first register. At the same time, executing pauli-X on the second register. Then the Hadamard transform on the second register is carried out, which results in $\frac{1}{2}\left({\left|0\right\rangle_{\rm{1}}{\rm{+}}\left|1\right\rangle_{1}}\right)\left({\left|0\right\rangle_{\rm{2}}{\rm{-}}\left|1\right\rangle_{2}}\right)\left({\left|0\right\rangle_{3}\left|{\varphi_{3}}\right\rangle_{45}+\left|1\right\rangle_{3}\left|{{{\vec{k}}}_{*}}\right\rangle_{4}\left|1\right\rangle_{5}}\right)$.

Step D3 Doing nothing when the first register is $\left|0\right\rangle_{1}$. Otherwise, performing swap operation for the second and the third register. Thus, the state can be obtained

Step D4 Measuring the first register with operator $M=\frac{1}{2}\left({\left|0\right\rangle-\left|1\right\rangle}\right)\left({\left\langle 0\right|-\left\langle 1\right|}\right)$, where the inner product of Eq. (14) and $\left|{{{\vec{k}}}_{*}}\right\rangle_{4}\left|1\right\rangle_{5}$ is $\sum\limits_{j=0}^{M-1}{\frac{{c^{\prime}\alpha_{j}^{2}}}{{\lambda_{j}+\sigma_{M}^{2}}}}$. Thus, the probability of measurement is $p=\frac{1}{2}\left({1+\sum\limits_{j=0}^{M-1}{\frac{{c^{\prime}\alpha_{j}^{2}}}{{\lambda_{j}+\sigma_{M}^{2}}}}}\right)$. Here, $c’$ is a constant. We can get specific value of $\sum\limits_{j=0}^{M-1}{\frac{{c^{\prime}\alpha_{j}^{2}}}{{\lambda_{j}+\sigma_{M}^{2}}}}$ from probability. Therefore, predicting covariance value in quantum context is achieved.

In classical gaussian process regression, covariance function is generally ${\mathop{\rm cov}}\left({f\left({\vec{x}_{p},\vec{x}_{q}}\right)}\right)=k\left({\vec{x}_{p},\vec{x}_{q}}\right)=\exp\left({-{\textstyle{1\over 2}}\left|{\vec{x}_{p}-\vec{x}_{q}}\right|^{2}}\right)$. Thus, the representation of covariance matrix and kernel function vector are respectively

The preparation of quantum initial states has always been difficult in quantum machine learning. And in this paper, it is also a challenge. Inspired by quantum radial basis network [9], coherent state [26, 27] and block coding techniques [29, 30] are utilized to achieve the preparation of initial states. Coherent states play a crucial role in quantum optics and mathematical physics. They are defined in the Fock states $\left\{{\left|0\right\rangle,\left|1\right\rangle,...}\right\}$, which is a basis of the infinitely dimensional Hilbert space ${\rm H}$. Let’s $a,a^{\dagger}$ respectively, be the annihilation and creation operators of the harmonic oscillator. Then

For any $n\geq 1$, it is easy to see that

Let $r\in R$ be a real number, its coherent state is defined by

It is a unit eigenvector of a corresponding to the eigenvalue $r$, that is $a\left|{\varphi_{r}}\right\rangle=r\left|{\varphi_{r}}\right\rangle$. From Eq. (19) and Eq. (20), we also have $\left|{\varphi_{r}}\right\rangle=e^{-{{r^{2}}\mathord{\left/{\vphantom{{r^{2}}2}}\right.\kern-1.2pt}2}}e^{ra^{\dagger}}\left|0\right\rangle=e^{r\left({a^{\dagger}-\frac{a}{2}}\right)}\left|0\right\rangle$, thus $\left|{\varphi_{r}}\right\rangle$ is obtained by a unitary operator of dimension infinity.

As for the preparation of $\left|{\varphi_{r}}\right\rangle$ in a finite quantum circuit, we can consider its Taylor approximation:

Thus, the upper bound on error square is

By stirling formula, we can get $T!\approx\sqrt{2\pi T}\left({\frac{T}{e}}\right)^{T}$. Now, keeping error within $\delta$. Let’s take Eq. (17) $\leqslant\delta^{2}$, that is $T$ only need to satisfy $2\log\frac{1}{\delta}-\frac{1}{2}\log 2\pi\leqslant\left({T+\frac{1}{2}}\right)\log T-2T\log r-T$. Thus, it is easy to get the Taylor approximation of $\left|{\varphi_{r}}\right\rangle$.

From the previous analysis, $\left|{\varphi_{r}}\right\rangle$ can be obtained by an unitary operator applying to $\left|0\right\rangle$. Because annihilation operators and creation operators are Hermitian, operator $e^{r(a^{\dagger}-\frac{1}{2}a)}$ can be simulated by LCU. It is also shown that $e^{r(a^{\dagger}-\frac{1}{2}a)}$ can simulate it with finite dimensions, which get a similar result. Thus for any vector ${{\vec{x}=}}\left({x_{1},x_{2},...,x_{N}}\right)^{T}$, the definition of coherent state is $\left|{\varphi_{x}}\right\rangle=\left|{\varphi_{x^{(1)}}}\right\rangle\otimes\left|{\varphi_{x^{(2)}}}\right\rangle\otimes\cdots\otimes\left|{\varphi_{x^{(N)}}}\right\rangle$, which costs $O\left({\frac{{Nt{\text{poly}}\log N\log\left({\frac{t}{\varepsilon}}\right)}}{{\log\log\left({\frac{t}{\varepsilon}}\right)\delta}}}\right)$. Now, considering the superposition states of training samples’ coherent states:

Firstly, applying the Fourier transform on $\left|0\right\rangle_{e1}$, which results in $\frac{1}{{\sqrt{M}}}\sum\limits_{p=0}^{M-1}{\left|p\right\rangle_{e1}}$. Then appending one qubit, which gets $\frac{1}{{\sqrt{M}}}\sum\limits_{p=0}^{M-1}{\left|p\right\rangle_{e1}\left|0\right\rangle_{e2}}$. Finally, applying unitary opeartors to the e2 register, Eq. (26) can be attained in time $O\left({\frac{{MNt{\text{poly}}\log N\log\left({\frac{t}{\varepsilon}}\right)}}{{\log\log\left({\frac{t}{\varepsilon}}\right)\delta}}}\right)$.

Then, taking the partial trace on the e2 register of $\left|\Psi\right\rangle\left\langle\Psi\right|$ gives rise to the density operator of covariance matrix

For covariance vector $\sum\limits_{p=0}^{M-1}{\exp\left({-\frac{1}{2}\left|{\vec{x}_{*}-\vec{x}_{p}}\right|^{2}}\right)}\left|p\right\rangle$, we put the target vector ${\vec{x}_{*}}$ into the training dataset D, which make dataset become $\left({\vec{x}_{i},y_{i}}\right)_{i=1}^{M+1}$. Here ${\vec{x}_{M+1}}$ is ${\vec{x}_{*}}$ and $y_{i}=0$. Then, computing new covariance matrix $K^{\prime}$ that $K^{\prime}=\frac{1}{{M+1}}\sum\limits_{p,q=0}^{M}{\exp\left({-\frac{1}{2}\left|{\vec{x}_{p}-\vec{x}_{q}}\right|^{2}}\right)}\left|p\right\rangle\left\langle q\right|$.

Ref. [29] presented a definition that supposed that $K^{\prime}$ is an $m$-qubit density operator, $\gamma,\varepsilon^{\prime}\in R^{\dagger}$ and $k\in N$. Then there is a $(m+k)$-qubit unitary $U^{\prime}$ is $\left({\gamma,k,\varepsilon^{\prime}}\right)$ -block-encoding of $A$ so that $\left\|{A-\gamma\left({\left\langle 0\right|^{\otimes k}\otimes I_{2m}}\right)U^{\prime}\left({\left|0\right\rangle^{\otimes k}\otimes I_{2m}}\right)}\right\|\leq\varepsilon^{\prime}$, where $\left\|\cdot\right\|$ is the spectral norm, thus, $U^{\prime}=\left({\begin{array}[]{*{20}c}{A^{\prime}}&\cdot\\ \cdot&\cdot\\ \end{array}}\right)$ satisfied $\left\|{A^{\prime}-\gamma A}\right\|\leq\varepsilon$. Therefore, $A^{\prime}\left|q\right\rangle$ can be obtained by considering $U^{\prime}\left|0\right\rangle\left|q\right\rangle$. That is $U^{\prime}\left|0\right\rangle\left|q\right\rangle\approx\left({\begin{array}[]{*{20}c}A^{\prime}&.\\ .&.\\ \end{array}}\right)\left({\begin{array}[]{*{20}c}{\left|q\right\rangle}\\ 0\\ \end{array}}\right)=\left({\begin{array}[]{*{20}c}{A^{\prime}\left|q\right\rangle}\\ .\\ \end{array}}\right)=A^{\prime}\left|q\right\rangle\left|0\right\rangle+\left|0\right\rangle^{\bot}$. Then, a projection measurement is performed in the basis $\left\{{\left|0\right\rangle,\left|1\right\rangle}\right\}$, the probability of obtaining measurement outcome $A\left|q\right\rangle$ is $tr\left({A^{\prime 2}}\right)$. Here $A$ is a density operator, thus $tr\left({A^{2}}\right)\leqslant 1$. If $A$ is a pure state, $tr\left({A^{2}}\right)=1$. Otherwise, $tr\left({A^{2}}\right)<1$. When the eigenvalues of $A$ are small, we put a coefficient in front of the matrix $A$ in order to obtain $A\left|q\right\rangle$ with a high probability. Here Gilyén et.al [31] proposed that the block-encoding of $A$ is $\left({G^{\dagger}\otimes I}\right)\left({I\otimes SWAP}\right)\left({G\otimes I}\right)$, where $G$ is an unitary operator that get a purification $\left|0\right\rangle\left|0\right\rangle\to\left|\Psi\right\rangle$ by $G$ appling on $\left|0\right\rangle\left|0\right\rangle$. Thus, $K^{\prime}\left|M\right\rangle=\frac{1}{{M+1}}\sum\limits_{p=0}^{M}{\exp\left({-\frac{1}{2}\left|{\vec{x}_{M}-\vec{x}_{p}}\right|^{2}}\right)}\left|p\right\rangle=\frac{1}{{M+1}}\sum\limits_{p=0}^{M}{\exp\left({-\frac{1}{2}\left|{\vec{x}_{*}-\vec{x}_{p}}\right|^{2}}\right)}\left|p\right\rangle$ can be obtained. In conclusion, covariance matrix $K$ and kernel function vector $\vec{k}_{*}$ can be prepared. Finally, the preparation of $\left|{\vec{k}_{*}}\right\rangle$ needs $O\left({\frac{{MNt{\text{poly}}\log N\log\left({\frac{t}{\varepsilon}}\right)}}{{\log\log\left({\frac{t}{\varepsilon}}\right)\delta}}+T}\right)$. $O\left(T\right)$ always refer to the complexity of implement the blocking-encoding. Without loss of generality, $O\left(T\right)=O\left({\log M}\right)$.