We address the problem of nonlinear dynamical system identification in state space formulation using an approximate Gaussian process (GP) regression framework where the basis for GP are learned from data. Approximate GP inference is used to address the high computational cost of exact GP frameworks, and is based on basis function expansion concept. We address the design of these basis functions using data and eigenvalue decomposition framework. The need for such design arises in practical applications where the collected data has outliers, discontinuities and other non-stationary characteristics which limit not only the applicability but also the performance of traditional GP framework with inherent stationary assumptions. Using the available data set, a kernel eigenvalue problem is formulated which is solved using Monte Carlo techniques to construct basis functions. Then, a low-rank approximation to the exact GP is obtained by an approach similar to kernel principle component analysis (k-PCA). The coefficients in the basis function expansion and other unknown parameters are learned from data using sequential Monte Carlo technique. The proposed GP framework for dynamical system identification is tested and validated against fixed basis framework using simulations and experimental data. INDEX TERMS Gaussian processes regression, non-stationary kernel, kernel PCA, sequential Monte Carlo, kernel eigenvalue problem.