Gaussian Processes (GPs) are a powerful tool for regression tasks, but their computational complexity is a limiting factor. In particular, it must calculate the Cholesky decomposition of a large covariance matrix, which results in a lower triangular matrix and is used for making predictions and calculating the uncertainty of those predictions. This step has a cubic time complexity of and is the prevalent part of the computation. Therefore, we use task-based parallelism with the High Performance ParalleX (HPX) library and tiled algorithms that view the covariance matrix as smaller tiles of equal size. The tiled algorithm executes parallel tasks that each perform a Basic Linear Algebra Subprograms (BLAS) operation on those tiles to speed up the computation. This led to the creation of the Gaussian Process Regression using Asynchronous Tasks (GPRat) Python library, which implements GPs based on this approach with the HPX runtime system. GPRat aims to provide a scalable and efficient implementation of GPs for large datasets with an implementation in C++ while leveraging the simplicity of Python for ease of use and seamless interaction with other Python libraries. It utilizes Intel oneAPI Math Kernel Library (oneMKL) to perform the BLAS operations on the Central Processing Unit (CPU). Graphics Processing Units (GPUs) are well-optimized for parallel processing, in particular for BLAS operations that operate on large matrices. This thesis extends the GPRat library with accelerator support for NVIDIA GPUs to further speed up the Cholesky decomposition and the other GP operations that utilize its result. The CUDA library includes the cuBLAS library and cuSolver library which provide the necessary Linear Algebra PACKage (LAPACK) and BLAS operations for executing computations on GPUs. This thesis presents the addition of CUDA support to the GPRat library. In particular, it describes the implementation of the tiled Cholesky decomposition on GPUs and also the other GP operations that utilize its result and therefore perform its additional operations on the GPU as well. We then compare the performance of the GPU implementation to the existing CPU implementation on a single node including a dual-socket AMD EPYC 9274F CPU and an NVIDIA A30 24 GB GPU. The results show that the GPU implementation is faster when working with at least 128 training samples and fewer but larger tiles. We also compared the performance of the GPRat library in computing the full covariance matrix prediction to the CPU implementation, but also to the GPyTorch and GPflow libraries on both CPU and GPU, which are two popular and highly efficient GP libraries for Python. We find that the GPRat library is the fastest on the CPU, competitive with the other libraries on the GPU, and scales well with the number of training samples. However, GPflow is slightly faster than GPRat on the GPU for higher numbers of training samples, and the optimal number of matrix tiles gives hardly any performance benefit on a single node.