LARGE-SCALE ESTIMATION OF VARIANCE AND COVARIANCE COMPONENTS

This paper concerns matrix computations within algorithms for variance and covariance component estimation. Hemmerle and Hartley [Technometrics, 15 (1973), pp. 819-831] showed how to compute the objective function and its derivatives fbr maximum likelihood estimation of variance components using matrices with dimensions of the order of the number of coefficients rather than that of the number of observations, Their approach was extended by Corbeil and Searle [Technometrics, 18 (1976), pp. 31-38] for restricted maximum likelihood estimation. A similar reduction in dimension is possible using expectation-maximization (EM) algorithms. In most cases, variance components are assumed to be strictly positive. We advocate the use of a modification that is numerically stable even if variance component estimates are small in magnitude. For problems in which the number of coefficients is large, Fellner [Proc. Statistical Computing Section, American Statistical Association, 1984, pp. 150-154], [Comm. Statist. Simulation Comput. B, 16 (1987), pp. 439-463] discusses the use of sparse matrix methods for positive definite systems in EM algorithms. We show how to compute the likelihood functions and their derivatives via sparse matrix methods for symmetric-indefinite systems, thus making solution of a much wider class of large-scale problems realizable. Results are formulated for the more general case of covariance components whenever possible.