Some genomic evaluation models require creation and inversion of a genomic relationship matrix (G). As the number of genotyped animals increases, G becomes larger and thus requires more time for inversion. A single-step genomic evaluation also requires inversion of the part of the pedigree relationship matrix for genotyped animals (A(22)). A strategy was developed to provide an approximation of the inverse of GG˜(-1) that may also be applied to the inverse of A(22)(A˜(22)(-1)) The algorithm proceeds by creation of an incomplete Cholesky factorization (T˜(-1)) of G(-1). For this purpose, a genomic relationship threshold determines whether 2 animals are closely related. For any animal, the sparsity pattern of the corresponding line in T˜(-1) will thus gather elements corresponding to all close relatives of that animal. Any line of T˜(-1) is filled in with resulting estimators of the least-squares regression of genomic relationships between close relatives on genomic relationship between the animal considered and those close relatives. The G˜(-1) was computed as the matrix product (T˜(-1))(')D(-1)T˜(-1), where D(-1) is a diagonal matrix. Then, T(-1)G(T(-1))(') resulted in a new matrix that is close to diagonal and also needs to be inverted. The inverse of that matrix was approximated with the same decomposition as for approximation of the inverse of G(G˜(-1)) and the procedure was repeated in successive rounds of recursion until a matrix was obtained that was close enough to diagonal to be inverted element by element. Two applications of the approximation algorithm were tested in a single-step genomic evaluation of US Holstein final score, and correlation coefficients between estimated breeding values based on either real or approximated G(-1) were compared. Approximations came closer to G(-1) as the number of recursion rounds increased. Approximations were even more accurate and expected to be faster for A(22). Timesaving strategies are needed to reduce the computing time required for the algorithm.