Today, Physics-informed machine learning (PIML) methods are one of the efective tools with high fexibility for solving inverse problems and operational equations. Among these methods, physics-informed learning model built upon Gaussian processes (PIGP) has a special place due to provide the posterior probabilistic distribution of their predictions in the context of Bayesian inference. In this method, the training phase to determine the optimal hyper parameters is equivalent to the optimization of a non-convex function called the likelihood function. Due to access the explicit form of the gradient, it is recommended to use conjugate gradient (CG) optimization algorithms. In addition, due to the necessity of computation of the determinant and inverse of the covariance matrix in each evaluation of the likelihood function, it is recommended to use CG methods in such a way that it can be completed in the minimum number of evaluations. In previous studies, only special form of CG method has been considered, which naturally will not have high efciency. In this paper, the efciency of the CG methods for optimization of the likelihood function in PIGP has been studied. The results of the numerical simulations show that the initial step length and search direction in CG methods have a signifcant efect on the number of evaluations of the likelihood function and consequently on the efciency of the PIGP. Also, according to the specifc characteristics of the objective function in this problem, in the traditional CG methods, normalizing the initial step length to avoid getting stuck in bad conditioned points and improving the search direction by using angle condition to guarantee global convergence have been proposed. The results of numerical simulations obtained from the investigation of seven diferent improved CG methods with diferent angles in angle condition (four angles) and diferent initial step lengths (three step lengths), show the signifcant efect of the proposed modifcati