A KNN Based Kalman Filter Gaussian Process Regression

1 Introduction

Over the past years, Bayesian parametric approaches have been well-developed for regression analysis [Bishop, 2006]. However, the data flexibility usually cannot be established by parametric models in all realistic details. Hence, Bayesian nonparametric methods have become popular due to the plausibly probabilistic framework. Gaussian Process (GP) [Rasmussen and Williams, 2006], a nonparametric model with an elegant probabilistic form, has been applied for regression problems such as robotics [Plagemann et al., 2007; Ko and Fox, 2008; Deisenroth et al., 2011; Wang and Chaibdraa, 2012], geophysics [Paciorek and Schervish, 2004], econometrics [L´azaro-Gredilla and Titsias, 2011] and so on. However, there exist several difficulties in GP due to the data properties in the real world [Plagemann, 2008].

First of all, the size of real data sets N is usually quite large, and it leads to the unsatisfactory computational load in GP (O(N 3 )). A number of sparse forms have been proposed for computation reduction. Concretely, by introducing a small pseudo-input training data [Snelson and Ghahramani, 2006] or spectral representation of GP [L´azaro-Gredilla et al., 2010], sparse GPs efficiently solve several large data sets in robotics. However, the approximation accuracy still needs to be further improved since sparse GPs may fall into the overfitting risk. Recently, a Kalman filter Gaussian process (KFGP)

[Reece and Roberts, 2010] has been developed by correlating the latent function values of different data subsets in the wellknown Kalman filter (KF) framework to reduce computation load. But KFGP still has to pay the price of accuracy due to the random training subset selection.

Secondly, many real data sets in geophysics or biology reflect the input-dependent smoothness (nonstationarity) where the correlation between function values f(x) and f(x′ ) does not only depend on x-x′ but it is also related to input locations x and x′ themselves [Rasmussen and Williams, 2006]. Usually, it is hard for the standard GP to correctly interpret the nonstationarity. A direct method is to define a nonstationary covariance function of GP [Paciorek and Schervish, 2004]. However, the eigenvalue decomposition for the covariance matrix of each input density is complicated since it requires a great number of parameters even in the low dimensional cases. Several extensions have been made by using kernel adaptation [Lang et al., 2007] and local smoothness [Plagemann et al., 2008]. Moreover, there also exists an infinite mixture model to address computation burden and nonstationarity together [Rasmussen and Ghahramani, 2002], but its hierarchical structure is difficult for membership determination in the gating network [Paciorek and Schervish, 2004]. A more practical approach is a local mixture of Gaussian processes [Urtasun and Darrell, 2008] that contains offline training hyperparameters and online inference using local GP experts. However, its accuracy depends on the number of local GP experts.

In this paper, we propose a novel KNN based Kalman filter Gaussian process (KNN-KFGP) method for various regression tasks. We firstly apply a test-input driven KNN search to construct a small training subset for each test, and then establish a state space model by correlating GP prior for different training subsets. Finally, a KF is applied to infer the latent function values in a predict-update manner. Our approach reduces computational load since the estimation for each test is based on a KNN constructed training subset. Meanwhile, it preserves the accuracy since KF produces a statistically optimal solution for linear and Gaussian systems [Kalman, 1960]. Moreover, KNN provides a locally correlated training subset for each test point, which is a suitable way to deal with the nonstationarity.

The rest of this paper is organized as follows: In Section 2, we briefly review the standard GP regression. In Section

3, we firstly introduce a construction procedure for the KNN based state space model, then apply KF to infer latent function values. In Section 4, we evaluate our proposed KNNKFGP by using several synthetic and real data sets. Finally, we conclude the paper in Section 5.

2 Gaussian Process Regression

The goal of the nonlinear regression is to learn a function f(x) in the noisy relationship between an input x ∈ Rd and its corresponding output y ∈ R

where the noise ϵ is typically assumed as a Gaussian distribution ϵ ∼ N(0, σ2 ). Since parametric approaches restrict the richness of the assumed function family when solving the regression problem, Bayesian nonparametric methods have recently become more attractive by giving a prior probability to every possible function.

Gaussian process (GP) is a popular nonparametric model due to the fact that it provides a flexible manner to represent a distribution over functions. Formally, it is a collection of random variables, any finite number of which have a joint Gaussian distribution [Rasmussen and Williams, 2006]. Similar to a Gaussian distribution parameterizing by a mean vector and covariance matrix, a GP is fully specified as f(x) ∼ GP (m(x), k(x, x′ )) with its mean function

and covariance function

In this paper, we follow a commonly used setting [Rasmussen and Williams, 2006]. The mean function in (2) is set to be zero and the covariance function in (3) is chosen as a stationary squared exponential covariance function k(x, x′ ) = σ2 f exp[-0.5(x - x′ )T L-2 (x - x′ )] where L = diag([l1 · · · ld]) is the length-scale matrix. For convenience, all the hyperparameters are denoted into a vector θ = [σ σf l1 · · · ld]T .

GP solves a regression problem in two steps. Firstly, learning the θ given the training set (X, y) = {(xi, yi)}Ni=1 . To do that, we apply a gradient based optimizer as presented in [Rasmussen and Williams, 2006] by maximizing the marginal likelihood p(y|X, θ). Secondly, inferring the predictive distribution of the function value at test inputs X⋆ = {xi⋆ }Mi=1 . In the context of GP, the probability p(y, f(X⋆)|X,X⋆, θ) is a jointly Gaussian distribution because any finite set of function values is jointly Gaussian distributed and the noise ϵ in (1) is independent identically distributed (iid) Gaussian noise. Hence, based on the conditional property of Gaussian distribution, the predictive distribution p(f(X⋆)|X, y,X⋆, θ) is also Gaussian [Rasmussen and Williams, 2006]:

where Kθ(X⋆, X) denotes an M × N covariance matrix in which each entry is computed by the covariance function k(x, x′ ) with the learned θ. The Kθ(X, X) and Kθ(X⋆, X⋆) are constructed in a similar way.

3 KNN Based Kalman Filter Gaussian Process Regression

In practice, GP suffers from an unsatisfactory computational burden (O(N 3 )) if the size of data set N is beyond a few thousand. By constructing a state space model with several small data collections, the recent KFGP [Reece and Roberts, 2010] provided an elegantly predict-update KF framework for computation reduction (O(MK3 ), where M is the number of KF iterations and K is the maximum size of all the training subsets). However, the blindly random training data selection for each test in KFGP largely deteriorates the estimation performance. Hence, we firstly need to explore an efficient mechanism of training subset selection for the test points in order to improve the accuracy of Bayesian inference in the state space model.

3.1 KNNConstructed State Space Model

Since the goal is to accurately estimate latent function values at the test inputs {xi⋆ }Mi=1 , we propose a test-input-driven KNN method to explore the strongly correlated local structure of each test point. Concretely, for xt⋆ , we find its Kt nearest training inputs Xt according to the Euclidian distance, then use Xt and the outputs yt of Xt as a small training set (Xt, yt) = {(xi, yi)}Kt i=1 for xt⋆ .

We now use KNN based multiple data collections to establish a state space model including the state and observation equation. Firstly, the state equation represents the Markovian transition of latent function values. Hence it is necessary to explore the relation between the (t - 1)th and tth data subset. We denoted Xc t = Xt ∪ xt⋆ and f c t = f(Xc t ). If f(x) follows a GP, the prior p(f c t , f c t-1 |Xc t-1 , Xc t , θ) is Gaussian distributed:

Then the conditional distribution of f c t is also Gaussian:

where

This conditional probability (6) actually reflects that the state equation is a linear evolution between f c t and f c t-1 with an additive Gaussian noise vf t ∼ N(0,Q(θ)):

Secondly, the observation equation, which is a relationship between hidden function values and outputs, is directly obtained from (1):

whereHt = [IKt 0] is an indexmatrix so that the observation yt is only related to its corresponding Xt by Htf c t = f(Xt). The covariance R(θ) = σ2 I of the noise vy t is originally from ϵ in (1). To sum up, our state space model is fully specified by (9-10). The whole construction is shown in Figure 1.

3.2 Bayesian Inference by Kalman Filter

After constructing the state space model in Subsection 3.1, we can now consider a regression task as a state estimation problem using recursive Bayesian filtering techniques. Given the learned θ, Equation (9) and (10) in the state space model are linear with additive Gaussian noises. The well-known Kalman Filter (KF) provides an optimal solution according to maximum-a-posteriori (MAP) [Kalman, 1960].

The detailed Bayesian inference procedure is described as follows. Suppose that the (t - 1)th posterior distribution p(f c t-1 |y1:t-1, X1:t-1, x1:t-1 ⋆ , θ) is a Gaussian distribution N (f c t-1|t-1 , P c t-1|t-1 ). The first step is to incorporate the transition distribution p(f c t |f c t-1 , Xt-1:t, xt-1:t ⋆ , θ) (directly from (6)) with the (t - 1)th posterior to make the prediction for the tth hidden function values:

We denote p(f c t |y1:t-1, X1:t, x1:t ⋆ , θ) as N(f c t|t-1 , P c t|t-1 ), then the prediction step is fully expressed as follows:

The second step is to update the predicted estimation with the tth observed output yt:

where the likelihood p(yt|f c t , θ, Xt, xt⋆ ) is directly from (10), the prediction p(f c t |y1:t-1, X1:t, x1:t ⋆ , θ) is from (11), and the

Algorithm 1 KNN-KFGP Regression

• 1: Training the hyperparameter vector θ using gradient optimization with an extra data set {(xi, yi)}Si=1
• 2: for t = 1 toM do
• 3: Applying KNN based data construction to find the training data (Xt, yt) = {(xi, yi)}Kt i=1 for xt⋆
• 4: Using θ to specify G(θ), Q(θ), R(θ) in (7-8) and (10) for state space model construction
• 5: Kalman Predict Step: Equation (12-13)
• 6: Kalman Update Step: Equation (16-18)
• 7: end for

marginal distribution is an analytical integral:

Finally we denote p(f c t |y1:t,X1:t, x1:t ⋆ , θ) = N(f c t|t , P c t|t ), then put (10), (11) and (15) into (14), we obtain the recursive form for the update step with the Kalman Gain Γt :

The whole algorithm is summarized in Algorithm 1. We notice that, if the standard GP with stationary covariance function (we call it as the stationary GP) can suitably learn the model of the data set, there is a tradeoff of accuracy vs computation load. The stationary GP establishes the whole correlation between the full training set and each test point, so the accuracy is high. However, its computational load is also high. Our KNN Kalman filter Gaussian process (KNNKFGP) is computationally efficient, compared with the stationary GP. The main computation reduction is the estimation step that is mainly governed by Kalman filter O(MK3 ) and KNN search O(MN) after various precomputations [Urtasun and Darrell, 2008]. The K is the maximum value of Kt, M is the number of test inputs and N is the total number of training points. But the accuracy of our KNN-KFGP is not as good as the stationary GP in this case due to the fact that our KNN-KFGP selects a small training subset for each test point. However, its accuracy will not decrease much. The reason is that the small training subset in our KNN-KFGP is strongly-correlated to the corresponding test point, and the estimation of f c t-1 can be used to improve the estimation of f c t in the predict-update Kalman filter framework.

Moreover, if the data set reflects the nonstationary phenomenon, the stationary GP cannot correctly model it. On the contrary, our KNN-KFGP allows each test point to find its strongly correlated small training subset, which is a suitable way to capture the local data property to deal with nonstationarity problems.

The computational efficiency is evaluated by Running Time (RT). In the SMSE and MNLP, the predictive mean μi⋆ for the ith test output is directly obtained from (17) and the predictive variance σi⋆ is the sum of the learned noise variance σ2 and the function variance from (18) .

4.1 The PendulumData

The first data set is a realistic simulation of a mechanical pendulum. There are 315 training and 315 test input-output pairs in this data set. Both the training and test sets are randomly ordered. In each pair, the input is a 9-dimensional pendulum parameter vector, and its corresponding output is an angular velocity of the pendulum.

In this experiment, we compare our KNN-KFGP with GP, Kalman Filter Gaussian Process (KFGP) [Reece and Roberts, 2010], Sparse Pseudo-input Gaussian Process (SPGP) [Snelson and Ghahramani, 2006] and Sparse Spectrum Gaussian Process (SSGP) [L´azaro-Gredilla et al., 2010]. We then report the prediction performance as a function of K in Figure 2, where K respectively represents the size of training set for each test point in KFGP and KNN-KFGP, the size of pseudoinput set for each test point in SPGP, the number of the basis function pairs for each test point in SSGP. Moreover, we learn the hyperparameters using gradient optimization with the full training set. For each K, we run all the methods three times and calculate the average SMSE, MNLP and RT (seconds).

Since the prediction for each test point in the stationary GP is made with the whole training set, we use its results as the

4 Experiments

In this section, we evaluate our KNN-KFGP on several synthetic and real data sets. We aims to achieve two objectives in these experiments. Firstly, if the data sets can be successfully modeled by the stationary GP, we will evaluate whether our KNN-KFGP is a computationally efficient and accurate approximation of the stationary GP. The pendulum and kin40k data sets 1 are used to validate this objective. Secondly, if the data sets cannot be suitably modeled by the stationary GP, we will test whether our KNN-KFGP is able to deal with the nonstationarity. The mobile robot perception and global land-surface precipitation data sets 2 are chosen for this objective. For all the data sets, the accuracy is evaluated by the test Standardized Mean Square Error (SMSE)

and the test Mean Negative Log Probability (MNLP)

1 The pendulum set and the kin-40k data set are available at: http://www.tsc.uc3m.es/ miguel/downloads.php

2 The global land-surface precipitation data set is available at: http://www.cgd.ucar.edu/cas/catalog/surface/precip/gpcc.html

base line. We found that it successfully captures the data flexibility with the best SMSE and MNLP. However, the tradeoff is the highest computation load. On the contrary, SSGP is more computationally efficient than stationary GP. The SMSE of SSGP is also good when K is small. But its MNLP becomes poorly deteriorated as K increases. The reason is that SSGP falls into the overfitting trap [L´azaro-Gredilla et al., 2010]. Compared to SSGP, SPGP is more acceptable for this data set. Moreover, the accuracy of KFGP is slightly better than SPGP, but the computation load is higher than SPGP.

To sum up, there is a tradeoff between computation load and accuracy for all the methods above. It is worth mentioning that this pendulum data is small, so all the methods are competitively fast. However, in terms of accuracy, our KNNKFGP is the most suitable choice among all the approximate GPs. In order to convincingly show our KNN-KFGP outperforms other methods in terms of the tradeoff between accuracy and computational efficiency, it is necessary to further evaluate all the methods in the following intensive data set.

4.2 The Kin-40k Data

The second data set is a realistic simulation of a robot arm. The randomly ordered training (10000 instances) and test (30000 instances) set are massive. In each instance, the input is a 8-dimensional parameter vector, and the output is scalar. In Figure 3, we show the prediction performance for all the methods (stationary GP, SPGP, SSGP, KFGP and our KNNKFGP) as a function of K. The notation of K is the same as the pendulum experiment. Because this data set is massive, we learn the hyperparameters using gradient optimization with a randomly selected training subset including 1000 points. For each K, we run all the methods three times and calculate the average SMSE, MNLP and RT(seconds).

The stationary GP makes the prediction for each test point with the whole training set, so its accuracy results (SMSE and MNLP) are the best. However, the computation burden is unsatisfactorily high. In contrast, KFGP is computationally efficient while its accuracy is bad. Furthermore, the computation load of SPGP, SSGP, our KNN-KFGP is similar and much better than the stationary GP. But the best accuracy among all the approximation methods is obtained by our KNN-KFGP. Hence, compared to other approximate GPs, our KNN-KFGP is a more accurate and efficient approach for large data sets.

4.3 Mobile Robot Perception

Environment perception is a fundamental task for mobile robot systems. It is essential for a mobile robot to correctly interpret the sensor information in spatially nonstationary environments in order to find its precise location. Here we consider a perception task to evaluate how well our KNN-KFGP can deal with the nonstationarity. Concretely, a mobile robot is located at (0, 0) with π/2 heading. To find an exit, it uses its range sensors for obstacle detection. We simulate 200 training points where the input is a bearing angle in [0, π] and the output is the corresponding distance measurement. The test inputs are from 0 to π with an interval π/180. The observations at different sides of the exit show the discontinuity of the measured data, which is a case of spatial nonstationarity.

In this experiment, we compare our KNN-KFGP to the stationary GP, KFGP and local GP (LGP) which is a related work using local GP experts to capture the nonstationarity [Urtasun and Darrell, 2008]. The notation of K in KFGP and KNN-KFGP is the same as the first two experiments. In LGP, K is the number of local GP experts. It is worth mentioning that, for our KNN-KFGP, the accuracy will not be monotonically increasing in the nonstationary cases when K is increasing. The underlying reason is that, if the K is too large, some disturbed, even wrongly correlated training points will be selected for each test point. That is why the stationary GP cannot correctly model the nonstationarity. Moreover, the random training data selection in KFGP leads to its low accuracy. In order to improve the accuracy, KFGP has to use more training data for each test. Hence, its performance is prone to stationary GP which is poor for nonstationary cases.

In this experiment, We learn the hyperparameters using gradient optimization with the full training set, then run all the methods three times from K = 1 to 5. We found that KNNKFGP shows the best performance when K = 2, KFGP and LGP represents the best performance when K = 5. Notice that, we choose a range that contains the best K since we mainly focus on finding the best K for our KNN-KFGP. Of course, the accuracy of LGP is better when K is larger, but the computation load is also higher. The results of all the methods are illustrated in Figure 4 and Table 1. In Figure 4(a), the standard GP mistakenly predicts that there is no exit even though the robot has detected the data on both sides of the exit. In Figure 4(b), the prediction of KFGP fails since the randomly selected training points in this case are disturbed. In Figure 4(c), LGP correctly finds the exit, but the uncertainty at the exit is high. Finally, in Figure 4(d) our KNN-KFGP successfully models the exit since the training subset for each test is strongly correlated to that test. From Table 1, it is also shown that LGP and our KNN-KFGP are suitable approaches for nonstationarity modeling, but our KNN-KFGP performs better with higher accuracy than LGP.

4.4 Global Land-Surface Precipitation

In practice, the environmental data sets usually exhibit the spatial nonstationarity. Hence we choose the global landsurface precipitation (January 2010) to evaluate whether our KNN-KFGP can handle the real-life nonstationarity. There are 3306 points in the training set. For each point, the input is the location and the output is the precipitation value. The test set is based on a regular grid with a spatial resolution of 2.5o × 2.5o latitude by longitude. As before, we compare our KNN-KFGP to the stationary GP, KFGP and LGP. The notation of K is the same as the robot perception experiment. We train the hyperparameters using gradient optimization with a randomly selected training subset including 1000 points, then run all the methods three times from K = 1 to 10. We found that KNN-KFGP shows the best performance when K = 5, KFGP and LGP represents the best performance when K = 10. The results for all the methods are illustrated in Figure 5 and Table 2. In Figure 5, the prediction of KFGP almost fails due to the random training set selection. In contrast, our KNN-KFGP flexibly learns the land-surface precipitation model. The underlying reason is that each test point is learned by using its strongly correlated training set and the predict-update KF mechanism preserves the accuracy. From Table 2, it is also shown that LGP and our KNN-KFGP are suitable approaches for this spatial nonstationary case, but our KNN-KFGP performs better than LGP. At last, it is worth mentioning that the RT of our KNN-KFGP is 475s which is much more efficient than GP (15374s).

5 Conclusions

In this paper, we have proposed a novel KNN-KFGP for regression. Firstly, small training subsets, which are constructed by the test-driven KNN search, are used to establish a GP prior based state space model. Then a KF is applied to infer the latent function values in a predict-update manner. Compared to other sparse GPs, our KNN-KFGP finds a more accurate solution with higher computational efficiency for large data sets. Moreover, the KNN search provides each test with its strongly correlated training set, so our KNNKFGP can more suitably solve the nonstationarity than GP. In the future, it will be interesting to explore an adaptive manner to learn K for each test to further improve the flexibility.