Quantum-Assisted Hilbert-Space Gaussian Process Regression

,


INTRODUCTION
Gaussian processes (GPs) are probabilistic machine learning methods widely used in applications such as robotics and control, signal processing, geophysics, climate modeling, financial markets, and data mining, as well as Bayesian optimization and probabilistic numerics [1][2][3][4].GPs are non-parametric probabilistic models that can be used for modeling multidimensional nonlinear functions through their mean and covariance functions [5,6].However, the traditional GP regression (GPR) methods struggle with computational efficiency, especially when handling large datasets [7].This limitation becomes particularly pronounced in fields where rapid processing of large-scale data is critical.In this paper, to tackle this challenge, we aim to accelerate the GPR through quantum computing.
The main computational complexity of GPR arises from the computation of the mean and variance of the posterior distribution.This process becomes increasingly computationally heavy with larger datasets, with computational and memory requirements scaling as O(N 3 ) and O(N 2 ), respectively, for N observations of input data.To alleviate this problem, various methods have been proposed.In the inducing point methods [8][9][10] the covariance matrix is approximated using a smaller number M of inducing points than the full training set, which reduces the computations to O(N M 2 ) or O(M 3 ) (for likelihood evaluation and prediction, respectively).
In this paper, we concentrate on low-rank methods [11][12][13] which are based on approximating the precision matrix via a set of M basis functions which also brings the computational complexity down to O(N M 2 ) or O(M 3 ).
In particular, we use the method proposed by Solin and Särkkä [7] which uses the Hilbert space of eigenfunctions defined by a Laplace operator to approximate the covariance function, which offers a tunable balance between computational complexity and approximation accuracy [7,14].
In recent years, quantum computers have emerged as potential replacements for classical computers [15].They offer exponential reductions in computational complexity for machine learning tasks.Quantum computing uses the principles of quantum mechanics to implement computational tasks and has been demonstrated for certain types of problems [16], for example, integer number factoring [17], fast database search [18], and matrix inversion [19].
Many quantum algorithms have been proposed for accelerating machine learning tasks.Among the plethora of quantum algorithms, the Harrow-Hassidim-Lloyd (HHL) matrix inversion algorithm [19] is often used to accelerate machine learning tasks.It serves as the foundation for various other algorithms such as quantum linear regression and quantum support vector machines [20][21][22].However, the HHL algorithm has its challenges, for example, in quantum state preparation, unitary simulation, sparsity, and matrix conditioning [23].
An HHL-based algorithm for quantum-assisted Gaussian-process regression was introduced in [24].This algorithm assumes that quantum state preparation and unitary simulation can be performed efficiently.Importantly, this algorithm addresses the inherent limitations of the HHL approach by appropriately selecting the covariance function to construct s-sparse matrices, and also carefully adjusts the noise parameters to ensure that the matrix remains well conditioned, as indicated by a condition number κ.To achieve a desired level of accuracy ϵ, it exhibits a run time that scales as O(log(N )κ 2 s 2 /ϵ).
The quantum principal component analysis (qPCA) is another quantum machine learning algorithm that draws inspiration from the HHL algorithm to estimate the dominant eigenvalues and eigenvectors [25].The authors in [26], proposed a method to prepare the covariance matrix on a quantum computer using annihilation and creation operators and implement the concept of qPCA to approximate the mean and variance of the GPR.This approach aims to achieve polynomial speedup compared to classical algorithms by overcoming the quantum state preparation and efficient unitary simulation assumptions.
The contribution of this paper is to integrate the Hilbert space approximation of the kernel presented in [7], into a quantum Gaussian process regression algorithm.This approach shifts the prediction complexity from O(N 3 ) to O(M 3 ), reducing the dependency from the number of observations N to the number of eigenfunctions M used to approximate the kernel.
Our methodology begins with the approximation of the kernel function using Hilbert space basis functions on a classical computer.Subsequently, we transfer this data matrix, characterized by a low-rank covariance function, onto a quantum computer.We then apply qPCA for extracting dominant eigenvectors and eigenvalues for non-sparse low-rank matrix into a quantum register.To derive the posterior mean and variance for reduced-rank Gaussian process regression, we employ conditional controlled rotations followed by the Hadamard tests for the mean and the Swap tests for the variance, respectively.We also include numerical examples to demonstrate and validate the effectiveness of our proposed method.Our proposed algorithm shows a polynomial speed advantage over existing classical algorithms for low-rank approximation in GP regression.
The structure of the paper is as follows.In Section II, we review the classical formulation for the Hilbert space approximation of GPR.We provide the quantumassisted Hilbert space GPR algorithm in Section III.The complexity analysis of the proposed algorithm and its comparison with state-of-the-art methods are given in Section IV.Section V discusses the numerical implementation of our algorithm on a classical simulator.We then conclude our findings in Section VI.

II. HILBERT SPACE APPROXIMATION OF GAUSSIAN PROCESS REGRESSION
In this section, we summarize the classical Hilbert space method for reduced-rank Gaussian process regression (GPR) [7] as well as show how GPR can be rewritten in terms of eigenvalues and eigenvectors.We first briefly review the classical GPR.Then, we show how to approximate the kernel using a Hilbert space of functions defined by the eigenspace of the Laplace operator.Finally, show how to express GPR in terms of singular value decomposition (SVD).This allows us to write these quantities in a suitable form so that they can be calculated using quantum states.

A. Gaussian process regression
Gaussian process regression [5] is a method for modeling and predicting multi-dimensional data.Consider a dataset D = (x i , y i ) N i=1 , where each {x i } N i=1 is a ddimensional input vector and y i is its corresponding mea-surement.In GPR, we aim to estimate an underlying function f (x) by modeling it as a realization of a Gaussian process.The measurements are then Gaussian distributed with added Gaussian noise ε i ∼ N 0, σ 2 : where k (x, x ′ ) denotes the covariance function (kernel), which is a positive semidefinite function k : Ω × Ω → R. The choice of kernel function drives the quality of the estimation.A common kernel choice for GPR is the square exponential covariance function [5]: where σ f and l are the signal scale and length scale hyperparameters respectively.
The objective in GPR is to predict the mean and variance of the output for new inputs x * .These predictions are derived from the posterior distribution, which is also Gaussian: The mean and variance of the posterior distribution are given by [5] Here, we denote by y the vector with components y i from the dataset, K the N × N matrix with entries K ij = k(x i , x j ) consisting of covariance functions between all input points in the training set, and k * is the covariance vector with the ith entry being k(x * , x i ).The kernel function can be approximated by a set of basis functions in a suitable Hilbert space as will be discussed next.

B. Kernel function approximation
We can approximate a kernel function by considering the eigenvalue problem of the Laplace operator [7]: where the domain Ω behaves well enough so that the eigenfunctions and eigenvalues exist.The functions ϕ j (•) are orthonormal with respect to the inner product which also defines a Hilbert space.All the eigenvalues λ j of the Laplace operator are real and positive.If the kernel function is isotropic k(x, x ′ ) = k(||x − x ′ ||) then its eigenvalues are given by the scalar function S(ω), called the spectral density, which is the Fourier transform of h → k(||h||).It turns out that we can approximate the kernel function in the domain Ω by [7] k Using this Hilbert space approximation of the kernel function, we can reformulate the Eqs.( 5) and ( 6).This modification allows for computationally efficient approximations for the mean and covariance of the GP: where Λ is a diagonal matrix with components Λ jj = S( λ j ), the matrix Φ has components Φ ij = ϕ j (x i ) and ϕ * has components ϕ j (x * ).We refer to this approximation as Hilbert space approximation for Gaussian process regression (HSGPR) [14].The approximation of the kernel now depends on the domain Ω and the set of eigenfunctions chosen in this domain.For the implementation of this paper, we chose Ω in the domain [−L, L].The Laplace operator in this domain gives rise to the set of sinusoidal eigenfunctions ϕ j (x) = L −1/2 sin(πj(x + L)/2L) with their corresponding eigenvalues λ j = (πj/2L) 2 .This kernel approximation allows us to reduce the complexity of the matrix inversion needed to find the mean and variance of the GPR.

C. Mean and variance of reduced rank GPR using singular value decomposition
In this section, we will convert the mean and variance expressions of GPR into a form that enables them to be expressed as expected values of quantum states and calculated in a quantum computer.Before applying our quantum algorithm, we modify Eqs.(10) and (11).For the GPR, we need the eigenvalues and eigenvectors of Φ ⊤ Φ + σ 2 Λ −1 which we wish to express in terms of Φ ⊤ Φ.We need to reformulate in such a way that both quantities share the same set of eigenvectors, allowing us to write the mean and variance of the GPR in terms of this common set of eigenvectors.This will enable us to write these quantities in terms of the expected values of quantum states.
To address this, we define where Now the eigenvectors of X ⊤ X + σ 2 I are the same as those of X ⊤ X.
We then begin by applying the SVD to the real data matrix X which is then expressed as X = UΣV ⊤ .Here Σ ∈ R R×R is a diagonal matrix containing the real singular values λ 1 , λ 2 , . . ., λ R and the orthogonal matrices U ∈ R N ×R (and V ∈ R R×M ) correspond to the left and right singular vectors, respectively.Taking into account the sum of X ⊤ X and σ 2 I, we derive Then, the eigendecomposition of X ⊤ X + σ 2 I −1 X ⊤ is given by: where i +σ 2 .We can write the Eq. ( 14) as Then, the mean of the GPR can be expressed using the SVD as Similarly, we can write the variance of GPR using the SVD as We have now expressed the mean and variance of GPR in a form that allows us to compute each of them as expected values of quantum states, which we will do in the next section.

III. QUANTUM-ASSISTED HILBERT SPACE GPR ALGORITHM
In this section, we propose a low-rank method for quantum-assisted Gaussian process regression which we call quantum-assisted Hilbert-space Gaussian process regression (QA-HSGPR).For its implementation, we have to encode the data matrix X ⊤ X into a quantum state.After that, we can implement a quantum algorithm that allows us to extract its eigenvalues.Then, we build the quantum circuits whose expected values correspond to the mean and variance that characterize the GPR.

A. Quantum state preparation from dataset
Quantum computers encode classical information into quantum states using qubits [27].A quantum state with FIG. 1.In this figure, qPCA is first employed on the matrix ρ X ⊤ X Following this, a conditionally controlled unitary operation is executed based on the eigenvalues register.Finally, we revert the additional τ qubit register to its original state by executing the corresponding inverse quantum operations to prepare the quantum state |ψ1⟩.
n qubits can be expressed as a 2 n dimensional vector |ψ⟩ = The coefficients a i are complex numbers that satisfy the normalization condition We use the notation ⟨ψ| to represent the conjugate transpose of the quantum state |ψ⟩.
We use an amplitude state encoding scheme to prepare the quantum state [28].The amplitude quantum state encoding encodes the classical vector (α into the coefficients of the quantum state.We begin by obtaining a matrix X ∈ R N ×M using the eigenfunction of the Laplace operator in the given domain.We then vectorize the matrix X and encode it using amplitude state encoding scheme [28] as Here, x m n represents the value of the classical data at position n, m of the data matrix X.It is important to note that the entries x m n must satisfy the condition n,m |x m n | 2 = 1.This ensures that the quantum state is properly normalized. Encoding data into a specific quantum state |ψ X ⟩ as given in Eq. ( 18) generally involves a computational complexity of O(N M ) in the conventional quantum state preparation methodologies [28,29].This complexity measure refers to the total number of quantum gates necessary to achieve the intended outcome.An approximate quantum amplitude encoding procedure has recently been proposed for more efficient state preparation [30].In that scheme, quantum state preparation is achieved in O(poly(log(N M ))), when dealing with a real data matrix.The low-rank kernel matrix X for HSGPR consists of real-valued entries, which would allow us to prepare it efficiently.

B. Estimation of eigenvalues
In this section, we show how to extract the eigenvalues of the symmetric matrix X ⊤ X and store them in an ancillary quantum register.This allows us to easily perform the conditional rotation operation, which is necessary to obtain the desired amplitude quantities in Eq. ( 16) and Eq.(17).Using the Gram-Schmidt decomposition of Eq. ( 18), we can reexpress |ψ X ⟩ as [31] |ψ Let us consider the density matrix ρ X ⊤ X = Tr n |ψ X ⟩ ⟨ψ X | by disregarding the |n⟩ register where Tr n is the partial trace on the n qubits, which can be written as Next, we apply the unitary evolution technique of quantum principal component analysis (qPCA) [25] ρ X ⊤ X to the register |m⟩ of |ψ X ⟩, resulting in for some large Z, and the states |ξ i ⟩ are intermediate states along the algorithm.By utilizing the quantum phase estimation algorithm, we can take the R dominant eigenvalues of the operator ρ X ⊤ X and write (cf.[31]) in which the singular values λ r are encoded in the τ qubits of an extra register.

C. Mean of Gaussian process regression
In this section, we provide the quantum method for computing the mean of GPR.We employ the conditional unitary on the ancilla qubit to invert the singular values.We add an extra ancilla qubit.The added ancilla qubit is conditionally rotated based on the eigenvalues register such that where the parameter c 1 is chosen such that the quantity c1 λ 2 r +σ 2 remains upper bounded by 1.After the conditional unitary, we reverse the computation in the τ qubits register by performing inverse operations of qPCA to bring them back into |0⟩ states The quantum circuit for preparing the quantum state |ψ 1 ⟩ is shown in Fig. 1. We Applying the Hadamard gate on the ancilla qubit leads to (24) Both |ψ 1 ⟩ and |ψ 2 ⟩ real vectors, and their inner products ⟨ψ 1 2 ⟩ and ⟨ψ 2 |ψ 1 ⟩ are equal.When measuring the ancilla qubit, the probability p(0) of measuring the ancilla in state 0 is given by where Thus we obtain an expression equal to the GPR mean as given in Eq. ( 16) up to a multiplicative constant.This mean value approximates the output function based on the data points and can be estimated using a quantum circuit.
In this figure, we illustrate the application of qPCA on the matrix ρ X ⊤ X .Initially, qPCA identifies the eigenvalues and eigenvectors of the matrix.We then apply a conditionally controlled unitary on ancilla register based on the eigenvalue register.Finally, the Swap test is employed to estimate the variance of the GPR.

D. Variance of Gaussian process regression
In this section, we build the quantum circuit which computes the variance of the Gaussian process regressor.For that purpose, conditional rotation is applied on the eigenvalues register such that where the parameter c 2 is chosen such that the quantity c2 λr √ λ 2 r +σ 2 remains upper bounded by 1.We proceed with the algorithm for measuring the ancilla qubit and consider only the measurements in the state |1⟩.Then, discarding the eigenvalue register, ancilla register, and right eigenvector register results in the final state where the probability of acceptance is given by We use the Swap test to obtain the variance of GPR.
This is the same expression as we derived in Eq. ( 17), up to a multiplicative constant.We then multiply with the noise variance σ 2 to obtain the variance of Gaussian process regressor.Fig. 3 shows the circuit implementation for computing the variance.

IV. COMPLEXITY ANALYSIS
In this section, we analyze the computational complexity associated with our proposed method.The algorithm starts with the quantum state preparation step.We employ an approximate quantum encoding scheme to prepare the quantum state |ψ X ⟩ ∈ R N ×M .This process requires a computational complexity of O(poly(log(N M ))).Similarly, the preparation of the quantum state |ψ 2 ⟩ mirrors this complexity.The total complexity for the preparation of the quantum state is O(poly(log(N M ))).
Following the state preparation step, we implement qPCA.The computational complexity for qPCA is O(log(M )ϵ −3 ) where ϵ denotes the desired error tolerance.The next phase involves a conditional unitary operation, achievable in O(log( 1 ϵ )).However, its complexity is relatively negligible compared with the complexity of qPCA.For both the mean and variance calculations in GPR, the initial algorithmic steps remain same.
To calculate the mean of the GPR, we employ the Hadamard test.The computational complexity of this test is linear in the number of qubits, with measurement accounting only for a constant factor which can be ignore.In the variance computation of the QA-HSGPR algorithm, the method involves measurement after the unitary conditional rotation.This requires O(κ 4 ) operations on average to measure the ancilla in the excited state.However, applying the techniques of [19,20], we can reduce this to O(κ 2 ).Following this, the Swap test is applied which is linear in the number of qubits.The measurement accounts for only a constant factor which can be ignored.Therefore, the overall computational complexity of the GPR is O(poly(log(N M )) log(M )ϵ −3 κ 2 ).A detailed comparison of the computational complexity of each step is summarized in Table I.
Classical Hilbert space methods for GPR generally endure a computational load of O(M 3 ) [7].In contrast, the mean and variance computations of our algorithm demonstrate a polynomially faster speed.We also compare our model with that of Zhao [24] whose algorithm complexity is O(log(N )ϵ −3 κ 2 ) dependent on the number of observations N assuming that the data matrix is already prepared in the quantum state.If we consider such an assumption, our method would have a complexity O(log(M )ϵ −3 κ 2 ), which is primarily dependent on the number of eigenfunctions M .This shifts the focus in complexity to M rather than N in our model, which significantly reduces the computational complexity, especially in scenarios with large datasets where usually M ≪ N .Furthermore, our method demonstrates significant improvements over the recently proposed quantum algorithm for Gaussian process regression.This contemporary model reports a time complexity of O(κ( 1 ) where P k denotes the probability of success for creating the quantum state and δ indicates the precision of the preparation of the state [26].This complexity depends on the dimension of the data points, which is not the case for our method.We present a detailed comparison of our method with existing approaches in Table II.This comparison reveals that the overall complexity of our proposed scheme is substantially lower than that of other existing methods.

V. NUMERICAL EXPERIMENTS
In this section, we present the numerical results of our proposed scheme.Our focus is on demonstrating the effectiveness of the method, by performing simulations on a classical computer.We use the built-in Qiskit function for quantum state encoding in our simulation [32].

A. Quantum circuit simulation
Several factors influence the performance of our method.A key aspect is the time parameter in the qPCA algorithm, which we use to estimate the eigenvalues λ 2 r in the quantum register, as shown in Eq. ( 22).This estimation is done using the unitary operator U = e −iρ X ⊤ X t , where we define the time parameter as t = 2π/δ R .Following the phase estimation bounds detailed in [33], we Zhao [24] -O(log N ϵ −3 κ 2 ) Chen [26] O( 1 √ can assert that δ R > λ 2 max , where λ 2 max represents the largest eigenvalue of operator ρ X ⊤ X .But, for a good approximation of eigenvalues, δ R should be slightly greater than λ 2 max .It could also happen that the qPCA estimation algorithm gives two different approximations to the same eigenvalue.To avoid this problem, we checked the similarity between the different excited states after the qPCA algorithm and discarded the states that likely represent the same eigenvalue in our simulation.The distinguishing of the eigenvalues becomes better when we increase the τ qubit register.
We then select the dominant R eigenvalues.The selection of the dominant eigenvalues R is a critical factor here.In our demonstration, the selection is made such that the lowest of the Rth eigenvalues exceeds 0.01.Specifically, the probability p of finding the desired state, as outlined in [31] is bounded by It is important to note that a significant decrease in the smallest eigenvalue will proportionally decrease the probability of measuring the desired state, needing a higher number of shots for an accurate estimation.We define the constants c 1 as λ 2 r + σ 2 and c 2 as λ r λ 2 r + σ 2 in our experiments.
To precisely mirror the classical results using a quantum computer, a substantial number of qubits and a high number of shots are required.Furthermore, optimizing the hyperparameters is crucial for effective implementation of the algorithm.Our algorithm is ideally suited for fault-tolerant quantum computers.

B. Simulation results
To demonstrate the functionality of our method, we have successfully implemented it on a much smaller scale.Our simulation involved N = 16 data samples derived from an oscillating function in a symmetric length interval L = 2π, with additive white Gaussian noise σ = 0.1, length scale l = 1, and signal variance σ f = 1.5.We use the 10 6 number of shots for this experiment and τ = 13 qubits for the eigenvalue register.We implement the approximation using a set of sinusoidal eigenfunctions in the domain Ω = [−L, L] to approximate the kernel.First, we chose M = 4 and performed estimations for R = 1, 2, 3, 4. We show the behavior of the mean for different R. The open-source implementation of our simulation is available in reference [34].
Fig. 4 compares the GPR using the exponential kernel, its Hilbert space approximation, and the reduced rank approximation implementing a quantum circuit proposed in this paper.We can see how with R = 4 the estimation already follows the tendency of the data.However, the estimation result is not exact.There are several reasons for this, first, we are using a limited amount of qubits in the precision of the eigenvalues, which reduces the precision of the mean estimation.Moreover, along the circuit, we have to implement multiple times controlled gates of the unitary operator e ιρ X ⊤ X t as well as controlled rotations of small angles, which lead to numerical errors in the simulations.
We also performed another simulation with M = 8 and the same number of data points with a different function.We performed an estimation with R = 4, as illustrated in Fig. 5.The additive white Gaussian noise σ = 0.1, length interval L = 2, length scale l = 1, and signal variance σ f = 0.5.We use τ = 16 qubits for the eigenvalue register and 10 6 shots.As can be observed in Fig. 5, the estimation of the mean and variance of GPR through a quantum computer gives a close approximation of HS-GPR.These simulations demonstrate the effectiveness of our algorithm and how it could be implemented when fault-tolerant quantum computers are available.

VI. CONCLUSION
In this paper, we have introduced a novel quantumassisted Gaussian process regression (GPR) algorithm leveraging a low-rank representation of the GP.Our algorithm addresses the high computational demand of GPR, showcasing how quantum computing can significantly enhance the scalability and efficiency of GPR models.A novel element of our contribution is the incorporation of the Hilbert space basis function approximation [7] into the quantum computing paradigm.This integration leads to significant improvements in computational efficiency, particularly in terms of reducing the computa- FIG. 5. Mean and variance of GPR using the squared exponential kernel (gray solid), the Hilbert space approximation of the kernel with M = 8 eigenfunctions (black dashed), and our QA-HSGPR (blue line) with data points N = 16 (red cross).Each point in the blue line represents a simulation.The shaded areas around each approximation line indicate the 95% confidence intervals, providing a visual representation of the uncertainty associated with each method.We can see that our proposed scheme approximates the HSGPR method well with R = 4.
tional complexity compared to classical algorithms.We also provide numerical examples within a quantum setting, which show that the method also works in practice.
As for future work, probabilistic numerics techniques [3] provide a means to obtain probabilistic approximations for numerical integrals.A Bayesian quadrature treats the integral as a Gaussian process [35,36].Being based on Gaussian process regression, Bayesian quadrature is faced with a significant computational challenge when evaluating the integral.The present methodology provides a promising method to evaluate large-scale integrals using Bayesian quadrature on a quantum computer.
The algorithm proposed here is suitable for faulttolerant quantum computers, which makes its implementation in NISQ devices a challenge for further work.The complexity of the circuit is mainly dominated by the qPCA and quantum phase estimation algorithm, then, alternative versions of these algorithms can be considered to reduce the complexity of the circuit.For the qPCA algorithm, a hybrid classical-quantum approach that implements a variational circuit can be considered to reduce the depth of the circuit [37,38].The previous proposal would reduce the depth of the circuit but increase the classical resources needed to execute the algorithm.On the other hand, it has been shown that iterative approaches of the quantum phase estimation algorithm reduce the complexity of the circuit needed for this task [39,40].The implementation of iterative versions of the quantum phase estimation algorithm would reduce the complexity of the circuit needed to implement our method, enabling the possibility of implementing it in quantum hardware.

4 FIG. 4 .
FIG.4.Mean of the GPR using the squared exponential kernel (gray), the Hilbert space approximation of the kernel with M = 4 eigenfunctions (black dashed line), and our reduced rank approximation using a quantum circuit (blue lines) with N = 16 data points (red cross).The blue lines range over R = 1, 2, 3, 4 showing how taking a larger rank increases the accuracy of the estimation.
then prepare another quantum state |ψ 2 ⟩ = |X * ⟩ |y⟩ |0⟩ |1⟩, where |X * ⟩ = l x l * |l⟩ and |y⟩ = l y l |l⟩ are normalized quantum states that encode the X * and y vectors respectively.We use the Hadamard test to estimate the inner product between these two states.The circuit diagram of the Hadamard test is shown in Fig. 2. The implementation of the Hadamard test begins with the application of a Hadamard gate on the ancilla qubit.Depending on the state of the ancillary qubit, different quantum states are generated: |ψ 1 ⟩ for the state |0⟩ and |ψ 2 ⟩ for the state |1⟩.This results in the composite quantum state

TABLE I .
The time complexity of each step in the proposed method.

TABLE II .
The time complexity of our proposed algorithm against existing quantum and classical counterparts.