  Gaussian process regression has proven very powerful in statistics, machine
learning and inverse problems. A crucial aspect of the success of this
methodology, in a wide range of applications to complex and real-world
problems, is hierarchical modeling and learning of hyperparameters. The purpose
of this paper is to study two paradigms of learning hierarchical parameters:
one is from the probabilistic Bayesian perspective, in particular, the
empirical Bayes approach that has been largely used in Bayesian statistics; the
other is from the deterministic and approximation theoretic view, and in
particular the kernel flow algorithm that was proposed recently in the machine
learning literature. Analysis of their consistency in the large data limit, as
well as explicit identification of their implicit bias in parameter learning,
are established in this paper for a Mat\'ern-like model on the torus. A
particular technical challenge we overcome is the learning of the regularity
parameter in the Mat\'ern-like field, for which consistency results have been
very scarce in the spatial statistics literature. Moreover, we conduct
extensive numerical experiments beyond the Mat\'ern-like model, comparing the
two algorithms further. These experiments demonstrate learning of other
hierarchical parameters, such as amplitude and lengthscale; they also
illustrate the setting of model misspecification in which the kernel flow
approach could show superior performance to the more traditional empirical
Bayes approach.
