In many applied sciences, the main aim is to learn the parameters in the operational equations which best fit the observed data. A framework for solving such problems is to employ Gaussian process (GP) emulators which are wellknown as nonparametric Bayesian machine learning techniques. GPs are among a class of methods known as kernel machines which can be used to approximate rather complex problems by tuning their hyperparameters. The maximum likelihood estimation (MLE) has widely been used to estimate the parameters of the operators and kernels. However, the MLE-based and Bayesian inference in the standard form are usually involved in setting up a covariance matrix which is generally ill-conditioned. As a result, constructing and inverting the covariance matrix using the standard form will become unstable to learn the parameters in the operational equations. In this paper, we propose a novel approach to tackle these computational complexities and also resolve the ill-conditioning problem by forming the covariance matrix using alternative bases via the Hilbert–Schmidt SVD (HS-SVD) approach. Applying this approach yields a novel matrix factorization of the block-structured covariance matrix which can be implemented stably by isolating the main source of the ill-conditioning. In contrast to standard matrix decompositions which start with a matrix and produce the resulting factors, the HS-SVD is constructed from the Hilbert–Schmidt eigenvalues and eigenvectors without the need to ever form the potentially ill-conditioned matrix. We also provide stable MLE and Bayesian inference to adaptively estimate hyperparameters, and the corresponding operators can then be efficiently predicted at some new points using the proposed HS-SVD bases. The efficiency and stability of the proposed HS-SVD method will be compared with the existing methods by several illustrations of the parametric linear equations, such as ordinary and partial differential equations, and integro-differentia
