The objective of this study was to compare the accuracy of genomic predictions in final score for young Holstein bulls calculated from single-step GBLUP models with the regular G⁻¹ and approximated G⁻¹ (G⁻¹_ap) matrices. The G⁻¹_ap was calculated with recursions on a small subset of animals. The regular G⁻¹ has a quadratic memory cost and cubic computational cost as the number of genotyped animals increases, whereas G⁻¹_ap has a linear cost for animals outside the subset. The predictor data set consisted of 77,066 genotyped animals, 9,009,998 pedigree animals, and 6,384,859 classified cows born in 2009 or earlier. For calculation of G⁻¹_ap, 9,406 high accuracy bulls or 16,828 high accuracy bulls and cows were used as the small subset. Genomic predictions (GEBV2009) were calculated for predicted bulls that had no classified daughters in 2009 but did in 2014. The validation data set contained phenotypes and pedigree recorded up to March 2014. Daughter yield deviations (DYD2014) were calculated for the predicted bulls with at least 30 daughters in 2014 (n = 2,948). Coefficient of determination (R²), calculated from a linear regression of DYD2014 on GEBV2009, was 0.44 with the regular G⁻¹ and 0.45 with G⁻¹_ap for either subset. Genomic predictions using all available 569,404 (570K) genotypes were also calculated with G⁻¹_ap on 9,406 bulls. The computation was performed using 16 CPU cores and 61 G bytes of working memory. Setting up G⁻¹_ap took 1.8 h, setting up matrices associated with A₂₂⁻¹ took 7 min, and iterations took 4.6 h, resulting in 6.4 total hours. In contrast, BLUP computations without genotypes took 1.6 h in total. The R² value from 570K genotypes was similar to the result from GEBV2009. Genetic predictions can be obtained with substantially less computational cost but without loss of reliability using G⁻¹_ap. The single-step GBLUP with G⁻¹_ap is applicable to very large genotyped populations.