I guess my first question would be, why do you want to generate Henderson's
equations?  There are much better ways to go about these calculations.  The
lme4 authors have just submitted a paper on the details of the calculations
used in fitting linear mixed models in lmer.  I believe Ben or Martin will
be making it available on arxiv.org soon.

Even if you want to do things this way there are several ways in which the
calculations could be expedited.  A calculation of the form

XtX <- t(X) %*% X

is more effectively expressed as

XX <- crossprod(X)

It is not clear why you want a ginv of XX.  A lot of symbolic computation
goes on in R to try to ensure that X is not rank deficient.  It is more
effective to take a Cholesky decomposition of XX (or a QR decomposition of
X) and solve using the Cholesky factor.

	[[alternative HTML version deleted]]

