Gaussian processes are the cornerstone of statistical analysis in many application ar- eas. Nevertheless, most of the applications are limited by their need to use the Cholesky factorization in the computation of the likelihood. In this work, we present a matrix-free approach for comput- ing the solution of the maximum likelihood problem involving Gaussian processes. The approach is based on a stochastic programming reformulation followed by sample average approximation applied to either the maximization problem or its optimality conditions. We provide statistical estimates of the approximate solution. The method is illustrated on several examples where the data is provided on a regular or irregular grid. In the latter case, the action of a covariance matrix on a vector is computed by means of fast multipole methods. For each of the examples presented, we demonstrate that the approach scales linearly with an increase in the number of sites.