for this approach to be efficient from the perspective of flops, one needs to store the Cholesky decompositions of the covariance matrices of $B_1Z,\ldots,B_q Z$, which would greatly increase the memory requirements of the algorithm. For dealing with truly massive data sets, our long-term plan is to combine the randomized approach studied here with a composite likelihood by using the randomized algorithms to compute the gradient of (\[composite\]), thus making it possible to consider $A_j$’s and $B_j$’s of larger rank than would be feasible if one had to do exact calculations.

Section \[sec2\] provides a bound on the efficiency of the estimating equations based on the approximate likelihood relative to the Fisher information matrix. The bound is in terms of the condition number of the true covariance matrix of the observations and shows that if the covariance matrix is well conditioned, $N$ does not need to be very large to obtain nearly optimal estimating equations. Section \[sec3\] shows how one can get improved estimating equations by choosing the $U_j$’s in (\[ascore\]) based on a design related to $2^n$ factorial designs. Section \[sec4\] describes details of the algorithms, including methods for solving the approximate score equations and the role of preconditioning. 