Technical note: Computing strategies in genome-wide selection

Genome-wide genetic evaluation might involve the computation of BLUP-like estimations, potentially including thousands of covariates (i.e., single-nucleotide polymorphism markers) for each record. This implies dense Henderson's mixed-model equations and considerable computing resources in time and storage, even for a few thousand records. Possible computing options include the type of storage and the solving algorithm. This work evaluated several computing options, including half- stored Cholesky decomposition, Gauss-Seidel, and 3 matrix-free strategies: Gauss-Seidel, Gauss-Seidel with residuals update, and preconditioned conjugate gradients. Matrix-free Gauss-Seidel with residuals update adjusts the residuals after computing the solution for each effect. This avoids adjusting the left-hand side of the equations by all other effects at every step of the algorithm and saves considerable computing time. Any Gauss-Seidel algorithm can easily be extended for variance component estimation by Markov chain-Monte Carlo. Let m and n be the number of records and markers, respectively. Computing time for Cholesky decomposition is proportional to n(3). Computing times per round are proportional to mn(2) in matrix-free Gauss-Seidel, to n(2) for half-stored Gauss-Seidel, and to n and m for the rest of the algorithms. Algorithms were tested on a real mouse data set, which included 1,928 records and 10,946 single-nucleotide polymorphism markers. Computing times were in the order of a few minutes for Gauss-Seidel with residuals update and preconditioned conjugate gradients, more than 1 h for half-stored Gauss-Seidel, 2 h for Cholesky decomposition, and 4 d for matrix- free Gauss-Seidel. Preconditioned conjugate gradients was the fastest. Gauss-Seidel with residuals update would be the method of choice for variance component estimation as well as solving.