Abstract The problem of estimating covariance matrices in balanced multivariate variance components models is discussed. As with univariate models, it is possible for the traditional estimators, based on differences of the mean square matrices, to produce estimates that are outside the parameter space. In fact, in many cases it is extremely likely that traditional estimates of the covariance matrices will not be nonnegative definite (nnd). In this article we develop an iterative estimation procedure, satisfying a least squares criterion, that is guaranteed to produce nnd estimates of the covariance matrices, discuss the speed of convergence, and provide an example to show how the estimates change.