To determine the optimal set of hyperparameters of a Gaussian process based on a large number of training data, both a linear system and a trace estimation problem must be solved. In this paper, we focus on establishing numerical methods for the case where the covariance matrix is given as the sum of possibly multiple Kronecker products, i.e., can be identified as a tensor. As such, we will represent this operator and the training data in the tensor train format. Based on the AMEn method and Krylov subspace methods, we derive an efficient scheme for computing the matrix functions required for evaluating the gradient and the objective function in hyperparameter optimization.