Restricted maximum likelihood (REML) is widely regarded as the preferred procedure for estimating variance components in animal breeding problems. The size of the coefficient matrix, however, often leads to computational difficulties and many simplified algorithms, including diagonalization, have been proposed. Diagonalization of the mixed model equations coefficient matrix, augmented by the numerator relationship matrix, in different steps using Householder’s transformations and the QL algorithm is proposed. The transformations need only be performed once. Very large data sets can be handled with ease and once the transformations have been performed, there is no practical limitations to the number of iterations that may be performed. A numerical example illustrating the procedure is supplied and a FORTRAN program, based on this approach, is available.