The construction of surrogate models, such as radial basis function (RBF) and Kriging-based surrogates, requires an invertible (square and full rank matrix) or pseudoinvertible (overdetermined) linear system to be solved. This study demonstrates that the method used to solve this linear system may result in up to five orders of magnitude difference in the accuracy of the constructed surrogate model using exactly the same information. Hence, this paper makes the canonic and important point toward reproducible science: the details of solving the linear system when constructing a surrogate model must be communicated. This point is clearly illustrated on a single function, namely the Styblinski-Tang test function by constructing over 200 RBF surrogate models from 128 Latin Hypercubed sampled points. The linear system in the construction of each surrogate model was solved using LU, QR, Cholesky, Singular-Value Decomposition, and the Moore-Penrose pseudoinverse. As we show, the decomposition method influences the utility of the surrogate model, which depends on the application, i.e., whether an accurate approximation of a surrogate is required or whether the ability to optimize the surrogate and capture the optimal design is pertinent. Evidently the selection of the optimal hyperparameters based on the cross validation error also significantly impacts the utility of the constructed surrogate. For our problem, it turns out that selecting the hyperparameters at the lowest cross validation error favors function approximation but adversely affects the ability to optimize the surrogate model. This is demonstrated by optimizing each constructed surrogate model from 16 fixed initial starting points and recording the optimal designs. For our problem, selecting the optimal hyperparameter that coincides with the lowest monotonically decreasing function value significantly improves the ability to optimize the surrogate for most solution strategies.