1. Introduction
Analyzing spatial data is an important task in a wide variety of f fields such as geology, ecology, climatology, soci- ology and urban planning. Gaussian process regression [1], or which is known as Kriging [2] in geostatistics, is a repre- sentative method for analyzing and interpolating spatial data. In Gaussian process regression, a covariance function plays a crucial role to define its behavior. With spatial data, kernels that solely depend on distance between two locations, e.g. Gaussian kernels, are usually used for covariance functions, since closely located points are assumed to have similar output values. However, some closely located locations can have different output values, and distant locations can have correlated output values. For example, weather at coastal area would be different from inland area in the same city, people's cultural behavior in two cities divided by a river would be dissimilar, and seismic activity would be related in distant areas when they share a fault line. In other words, the observed input space would be different from its intrinsic space that reflects relationships among locations. Revealing the intrinsic space is beneficial to understand the given spatial data as well as to improve interpolation performance. In this paper, we propose the input-regularized Gaussian process latent variable model (IGPLVM) that visualizes the latent intrinsic space. The proposed model assumes that each observed input location has its own latent location, and the covariance function is determined by distance between the latent locations. Since neighbor relationships in the latent space should be similar to those in the observed input space as long as the output covariance of the given data is captured by the latent locations, latent locations are regularized to preserve the neighbor relationships. When the regularization is omitted, the proposed model becomes the Gaussian process latent variable model (GPLVM) [3], [4], which is an unsupervised method that learns latent locations using output values without input information. When the regularization is strong enough to perfectly preserve the neighbor relationships, it corresponds to Gaussian process regression, where input locations are assumed to be noise free. The proposed model adaptively uses input information so that output variables are modeled properly. The proposed model improves interpolation performance by flexibly defining the covariance function by adjusting latent locations. Since the proposed model learns a common latent intrinsic space shared by multiple output variables, it can be used for multi-task learning. The latent space learned by using output values in a task helps to interpolate output values in another related task, even if data are too sparse to learn latent locations with a single task. Note that, although we motivate the proposed model for analyzing spatial data, the proposed model is applicable to other data for any regression problems where input and output values are contained, such as time-series, spatio-temporal, high dimensional input, and multiple output data. The remainder of this paper is organized as follows. In Section 2, we present our task and introduce Gaussian process regression, on which the proposed model is based. In Section 3, we formulate the proposed model and provide its learning procedures. We review related work in Section 4. In Section 5, the effectiveness of the proposed method is demonstrated by experiments with real-world spatial data sets in terms of interpolation performance of multiple out- put values. Finally, we present concluding remarks and a discussion of future work in Section 6.
2. Preliminaries
Suppose that we have a set of input and output instances, {{xn, Yn }} N=1, where xn ∈ RL and yn ∈ RD. For example, in a case of spatial climate data, x n is th n th location vector such as latitude and longitude, and y n is an observed weather condition vector at the location such as temperature, humidity and precipitation. With standard regression methods, an output y n is as- sumed to be generated by mapping input x n using nonlinear functions,
where f d ( · ) is the nonlinear function for the d th feature, and ε ～ N (0,β−1) is a Gaussian noise. Gaussian process regression uses a Gaussian process for a prior distribution of the nonlinear function,
where m(x) is a mean function, which is often set to zero, and kernel function k(x,x') specifies the covariance of outputs between two locations as follows,
Since the kernel determines the behavior of the nonlinear functions, it is important to use an appropriate kernel for the given data. For real-valued input data, such as time-series and spatial data, kernels that are negatively correlated to distance between two locations || x – x' || are often used, such as Gaussian kernels. However, output values can be different even if two locations are close together, and they can be similar even if two locations are far apart.
3. Input-regularized Gaussian Process
Latent Variable Models
We propose the input-regularized Gaussian process latent variable model (IGPLVM), which is a method to obtain kernels that appropriately capture the output covariance between inputs by distorting the input space in the Gaussian process regression framework. The distorted input space reveals intrinsic characteristics of the input space. The proposed model assumes that an output y n is gener- ated from its latent intrinsic location z n ∈ R K , instead of input location x n , as follows,
The dimensionality of the latent space K can be different from that of the input space L . Then, the probability of the output observations Y = ( Y 1 , . . , Y N ) T given the latent locations Z = ( z 1 , . . , z N ) , by integrating out the
where K is the N × N covariance matrix defined by the kernel function k(z n , z n ′ ), and θ is the kernel hyperparam- eter vector. In this paper, we use a Gaussian kernel with an additive noise term,
where δ n m  =  1 if n  = m  and δ n m  =  0  otherwise, and θ = (α,ℓ,β) are the kernel parameters. Note that the GPLVM [3] finds latent locations Z that minimizes the negative log likelihood of (5),
where input information X = ( x 1 ,  ·  · , x N ) T  is not used. We assume that neighbor relationships in the latent space should be similar to those in the input space as long as the output covariance that is specified by the latent space is appropriate to the given data. The proposed method models the neighbor relationships by defining a probability of being selected as neighbors as defined in stochastic neighbor embedding [5], [6]. In the latent space, the probability that n selects n' as its neighbors is given by
where locations with small Euclidean distance ‖ z n − z n ′ in the latent space are likely to be selected as its neighbors. Similarly, in the input space, the neighborhood probability is given by
probabilities are matched. This is achieved by minimizing the following sum of Kullback-Leibuler divergences between the probabilities,
The proposed method finds latent locations that properly capture the output covariance while preserving the neighbor relationships by minimizing the following sum of (7) and (10),
where λ > 0 is a hyperparameter that controls how the neighbor relationships are preserved. When λ = 0, it corresponds to GPLVM. When λ = ∞ and dimensionality of the latent space is the same with the input space K = L , it corresponds to Gaussian process regression since the latent locations become the same with the input locations Z = X. The proposed method can be seen as a multi-task
learning method based on Gaussian process regression. Since the latent locations are learned by using all of the output variables, the learned covariance matrix can improve multi- task regression performance when the output variables are related. A local optimum solution of latent locations Z and kernel hyperparameters θ is obtained by minimizing (11) using gradient-based optimization methods such as the quasi- Newton method [7]. The gradients of the GPLVM term E Y with respect to a latent location are calculated by
using the chain rule. The gradients of the regularization term Ex with respect to a latent location are calculated by
where a sum of forces pulling or pushing z depends on difference of the neighborhood probability p(n'|n,X) – p(n'|n, Z). We can select the hyperparameter value λ by cross- validation on an interpolation problem of the output matrix Y. The cross-validation procedure is as follows. The ele- ments of the output matrix are randomly split into multiple subsets. While the elements in a subset are supposed to be missing, the latent locations and kernel hyperparameters are estimated with a fixed hyperparameter value λ. The interpolation performance of the learned model is evaluated by the rooted mean squared error for the missing elements. The error is averaged over different subsets, and a hyper- parameter value that achieved the lowest error is selected. When Ynd is missing, its estimated value is given by Ynd = k n K-1yd, where kn = (k(zn,z1),...,k(zn,zN)) T and yd = (yid, . . ,yNd) T, and they are calculated using the observed data. Note that even when all of the output vari- ables are missing with an instance, the latent location can be estimated by using its neighbor relationships. The proposed model can be seen as a probabilistic gen- erative model. An output value y nd is generated from a Gaussian distribution with mean f d ( z n ) using the latent location z n , where the nonlinear function f d ( · ) is generated from a Gaussian process prior. The input locations x n are not directly generated, but we can consider that the neigh- bors of the input locations are generated by a multinomial distribution, where the multinomial parameters are defined by latent locations Z as in (8). Although we considered that the input variables are real- valued, the proposed method is applicable to other types of input variables that can calculate similarity between two input locations to define neighbor relationships in (9). For example, we can use normalized tree kernels and graph kernels for tree and graph data, respectively, instead of normalized Gaussian kernels in (9).
4.  Related Work
With Gaussian process regression, the covariance function is usually defined parametrically, e.g. Gaussian kernels, and its hyperparameters are estimated, e.g. length-scale and noise variance in Gaussian kernels, by maximizing the marginal likelihood. However, the hyperparameter tuning would not be enough to capture the covariance structure when the intrinsic space is different from the observed input space. The proposed model can define the covariance function by flexibly tuning latent locations. A complex form of covariance functions is learned by using multiple kernels [8], by concatenation different kernels [9], and by modeling a spectral density with Gaussian mixtures [10]. These methods are not designed for revealing the latent intrinsic space. The manifold Gaussian process [11] is related to our work. It transforms input locations by multi-layer neural networks so as to model complex nonsmooth functions. Since the proposed model optimizes a latent location for each input location nonparametrically so as to preserve the neighbor relationships, it can be used for visualizing the latent intrinsic space. In addition, by tuning a single hyperparameter, the proposed model can become the GPLVM and Gaussian process regression, which are successfully used in a wide variety of applications  [12], [13], [14], [15]. The warped Gaussian process [16] transforms the output space. Gaussian process regression with input noise [17] assumes that each input location is independently corrupted by Gaussian noise, and does not give latent locations explicitly. On the other hand, with the proposed model, distortion of each location is not independent since latent locations are estimated so as to preserve the neighbor relationships. This enables us to estimate latent locations even without output values. The GPLVM [3] is an unsupervised method for imputating output values by assuming latent locations. The GPLVM is not designed for regression tasks, and it does use input information. The discriminative GPLVM [18] is a type of regularized GPLVMs. It is used for classification, and out- puts are required to be categorical variables. The Gaussian process latent variable set model [19] is a regression method that learns latent locations to model a flexible covariance function. It assumes bag-of-features input, and each feature has a latent location. On the other hand, the proposed model can take any input type, such as real-valued, tree, graph, bag- of-features data, as long as similarity can be calculated, and each instance has a latent location. Many visualization methods have been proposed, such as multidimensional scaling, principle component analy- sis, GPLVM, Isomap [20], Laplacian eigenmap [21] and stochastic neighbor embedding (SNE) [5]. They embeds high dimensional unlabeled data into a low dimensional space. Supervised visualization methods have been also proposed, such as canonical correlation analysis (CCA), kernel CCA [22], Fisher linear discriminant analysis and discriminant GPLVM. For example, CCA embeds input and output pairs into a low dimensional space so as to maximize the correlation. The proposed model is a supervised visual- ization method that reveals an intrinsic input space for better output interpolation. The proposed model can be seen as a combination of GPLVM and SNE, where GPLVM is used for modeling output values, and SNE is used for preserving the neighbor relationships of input locations. In particular, the objective function is simply defined by a weighed sum of the objective functions of GPLVM and SNE.
5. Experiments
5.1 Data
We evaluated the proposed method by using a real-world spatial data set: the comprehensive climate data of North America (NA) 1 . The NA data set consists of monthly climate reports from 1990 to 2002 [23], [24]. We used 16 output variables, such as carbon dioxide and temperature, which were interpolated on 2.5 × 2.5 degree grid with 125 locations. We normalized all output values to mean zero and unit variance, and conducted experiments with data for each month, where the input variables were latitude and longitude.
5.2 Comparing Methods
We compared the proposed IGPLVM (Proposed) with the following five methods: GP, MTGP, GPLVM, MF and KNN. GP is a Gaussian process regression method, which assumes that output values are determined by the input locations using nonlinear functions with Gaussian process priors. MTGP is a multi-task Gaussian process regression method [25], which learns relationships between output variables. We used the code provided by the authors 2 . GPLVM is a Gaussian process based nonlinear matrix imputation method [3]. MF is a matrix factorization method [26]. It imputates missing values by the product of two low-rank matrices. Both GPLVM and MF do not use input information. KNN is a k -nearest neighbor regression method, which estimates a missing value by the average value of its four neighbors. The dimensionality of the latent space with the proposed method, GPLVM and MF was set at K = 2, which is the same with that of the input space. With the proposed method, we selected a hyperparameter λ for each output variable by five-fold cross-validation from {0, 1, 10, 10 2 , 10 3 , 10 4 , 10 5 }. The latent locations were initialized by the input locations.
5.3 Results
We measured the effectiveness of the proposed method by interpolation tasks. Ten percent of output values were randomly selected as test data. The performance was evalu- ated by rooted mean squared error (RMSE). Table 1 shows the RMSE with the NA data set. The proposed method achieved the lowest average RMSE. This result indicates that the output covariance is properly modeled by distorting the input space with the proposed method. GP achieved low RMSE with output variables whose covariance is determined by the input locations, such as WET and DTR. GPLVM achieved low RMSE with output variables which can be esti- mated easily from other output variables. Since the proposed method can become the GP and GPLVM by controlling the hyperparameter for each output variable, the RMSE of the proposed method was low for all output variables. The computational time of the proposed method was 18 seconds with a one-month data using a computer with Xeon 7350 2.93GHz CPU. Figure 1 shows the original input space and the estimated latent space by the proposed method. When λ = 0, the latent locations were different from the input locations, since the latent locations did not regularized by the input locations at all. When λ = 10 5 , the latent space was the same with the input space, since the effect of preserving neighbor relationships became dominant. With λ = 10 and λ = 10 3 , the some neighbor relationships were preserved but some latent locations were transformed so as to model the output covariance. The latent locations of west coast were separated from the other locations (b,c,d), because west coast exhibits different weather from the other area.
6. Conclusion
In this paper, we have proposed a probabilistic model for discovering a latent intrinsic space. The proposed method is based on Gaussian processes, where a latent location is associated with each input location, and output values are determined by the latent locations. The latent locations are estimated so as to preserve the neighbor relationships as well as to capture the output covariance of the given data. Although our results have been encouraging, our framework can be further improved upon in a number of ways. Firstly, we would like to estimate the hyperparameter for controlling input regularization by using the variational Bayesian framework [27]. Secondly, we plan to estimate re- lationships among output variables using multi-task learning techniques [25]. With the proposed method a single latent space is shared by all output variables. By using estimated task relationships, we can obtain multiple latent spaces that capture the characteristics of individual output variables, and it leads to better interpolation performance for missing values. Finally, we would like to extend our method more scalable by using scalable Gaussian process techniques [28], [29].

Acknowledgments
A part of the research results has been achieved by "Research and Development on Fundamental and Utiliza- tion Technologies for Social Big Data," the Commissioned Research of National Institute of Information and Commu- nications Technology (NICT), JAPAN.