Tensor-valued data are being encountered increasingly more commonly, in the biological, natural as well as the social sciences. The learning of the unknown model parameter vector given such data,involves covariance modelling of such data, but this can be difficult owing to the high-dimensional nature of the data, where the numerical challenge in such modelling can only be compounded by largeness of the data set. Assuming such data to be modelled using a correspondingly high-dimensional Gaussian Process, the joint density of a finite set of such data sets is then a tensor normal distribution,with density parametrised by a mean tensor and the k covariance matrices. When aiming to model the covariance structure of the data, we need to estimate the covariance matrices. We present a new method in which we perform such covariance modelling by first expressing the probability density of the available data sets as tensor-normal, i.e. the likelihood (of the unknown matrix and tensor-variate parameters of the GP that the data is modelled using, given the data) is tensor-normal. We then invoke appropriate (vague) priors on these unknown parameters and express the posterior of the unknowns given the data. We sample from this posterior using an appropriate variant of Metropolis Hastings. In order to reduce computational burden, the mean tensor is estimated by maximum likelihood estimation in a pre-processing step. We perform empirical illustration of the method using large three-dimensional astronomical and economic data sets of size N times M times P. To begin with, we choose to work with the squared exponential covariance function, leading us to learn (from the data), the correlation length-scales. Although the difficulty of learning the covariance model is reduced by the undertaken steps, inference with MCMC on these large, high-dimensional data sets is still time and resource intensive. This motivates us to use an efficient variant of the Metropolis-Hastings algorithm--Transformation based MCMC--employed to perform efficient sampling from a high-dimensional state space. In further applications, other kinds of covariance functions will be discussed as well. Once we perform the covariance modelling of such a data set, we will learn the unknown model parameter vector at which a measured (or test) data set has been obtained, given the already modelled data (training data), augmented by the test data.