This dissertation delves into the advanced realms of Gaussian Processes (GPs) and Bayesian Optimization (BO), presenting novel methodologies that enhance their performance and applicability. GPs, as a principled probabilistic approach, are powerful in modeling complex and noisy functions due to their non-parametric nature and capability for uncertainty quantification. However, exact GPs become intractable for large datasets since the computational cost scales cubically with the size of the dataset. In particular, this dissertation focuses on improving variational GPs, which is able to handle large-scale data by sparsifying the model via inducing points and approximating the posterior. Despite advances, variational GPs still may require many inducing points (and significant computational costs) to achieve good accuracy, a gap this dissertation aims to bridge.This dissertation also studies efficient computational methods for Bayesian transformed GPs (BTG), which is particularly useful when the Gaussian assumption is not satisfied and data is limited. Furthermore, the dissertation explores BO as a method for optimizing complex and expensive objective functions, with an emphasis on its application in hyperparameter tuning. By leveraging the probabilistic modeling strengths of GPs, BO can efficiently traverse the hyperparameter space, thus reducing the need for extensive model evaluations. Through the introduction of novel algorithms and methodologies, this research not only enhances the performance of BTG and variational GPs but also broadens the scope of BO in hyperparameter tuning.