Monte Carlo simulations in geostatistics often require the generation of realizations of multi-dimensional Gaussian processes with prescribed mean and covariance over a specified grid Ω of sampling points in R:d In principle, given a factorization C = AAT of the process’ covariance matrix C, realizations can always be constructed from the product Aϵ where ϵ is a white noise vector with unit variance. However, this approach generally has a high computational cost when compared with techniques such as the spectral or turning bands methods. Nevertheless, matrix factorization has advantages in that it does not require any special structure, such as isotropy or stationarity of the covariance model, o r regularity of the sampling grid Ω. Therefore, in an effort to speed up the approach at the cost of some loss in accuracy, Davis (1987b) proposed approximating the square root A by a low order matrix polynomial in C. This paper explores the construction of such polynomial approximations in more detail

The paper first presents a systematic error analysis which justifies choosing the polynomial to be a best fit to the square root function on the unit interval in either the minimax metric or a particular weighted least squares metric. Given the well known difficulty of constructing arbitrary order minimax polynomial approximations, the paper suggests approximation by Chebychev polynomial expansions to give reasonable minimax solutions, plus approximation by quasi-Jacobi polynomial expansions to give optimal weighted least squares solutions. Rates of convergence are given for both expansions: for continuous covariance models these are restricted by the singularity of the square root function at the origin, but they improve dramatically if a nugget effect is included in the model. Thus the method may well be practical for generation of simulations in irregular problems, especially if high fidelity to the specified covariance model is not required.