When a Gaussian Markov random field (GMRF) is used as a metamodel of an unknown response surface for a discrete optimization via simulation (DOvS) problem, the hyperparameters of the GMRF are estimated based on a few initial design points in a large feasible solution space. Although the maximum likelihood estimators (MLEs) are most commonly adopted to estimate these hyperparameters, its computation time increases polynomially in the size of the feasible solution space. We introduce new generalized method of moments (GMM) estimators of the hyperparameters of GMRFs and their initial sampling schemes, and show they are consistent under some conditions. Unlike MLEs, the computation time for these GMM estimators does not depend on the size of the feasible solution space. We show empirically that the GMM estimators have smaller biases and standard errors than MLE for a wide range of initial simulation budget while requiring orders of magnitude smaller computation time.