Avoiding the explicit inverse of the pedigree relationship matrix among genotyped animals in approximate reliabilities of single-step genomic predictions

Technical Summary: Avoiding the Explicit Inverse of the Pedigree Relationship Matrix in ssGBLUP Reliabilities

Problem Statement
In single-step genomic BLUP (ssGBLUP), computing approximate reliabilities for estimated breeding values (EBVs) requires the inverse of the pedigree relationship matrix restricted to genotyped animals ($A^{-1}_{22}$). While the genomic relationship matrix ($G$) bottleneck has been addressed by the Algorithm for Proven and Young (APY), $A^{-1}_{22}$ remains a significant computational hurdle. Explicitly forming and inverting $A^{-1}_{22}$ demands quadratic storage and cubic computation time. As the number of genotyped animals exceeds a few tens of thousands, this approach becomes prohibitive, preventing the calculation of approximate reliabilities in large-scale genomic evaluations. Existing methods for solving ssGBLUP equations (e.g., using the Schur complement or Colleau's algorithm) have not been fully adapted or compared specifically for the context of approximate reliability computation.

Methodology
The study evaluated five strategies designed to avoid the explicit formation of $A^{-1}_{22}$ when computing approximate reliabilities. These methods were tested against a "dense inverse" benchmark (DenseInv) using data from the PROMEBO Angus breeding program. The dataset included 433,583 pedigreed animals, of which 24,086 were genotyped, across 14 traits in eight multi-trait analyses.

The five strategies are:

1. SchurPCG-P: Uses the Schur complement identity ($A^{-1}_{22} = A_{22} - A_{21}A^{-1}_{11}A_{12}$) to express the matrix implicitly through sparse blocks of the full pedigree inverse ($A^{-1}$). It employs the Preconditioned Conjugate Gradient (PCG) method to solve $n_2$ linear systems iteratively, utilizing a parallel sparse factorization of the non-genotyped block ($A_{11}$) via PARDISO.
2. WoodColl (Woodbury–Colleau): Reformulates the problem using the Woodbury identity to replace $A^{-1}_{22}$ with $A_{22}$ itself. It computes the diagonal of the resulting matrix by solving $n_2$ linear systems using the Colleau algorithm, which evaluates matrix-vector products ($A_{22}v$) without storing the dense matrix.
3. HutchColl (Hutchinson–Colleau): Combines the Woodbury reformulation with a stochastic Hutchinson estimator. Instead of solving $n_2$ systems, it uses a small number of random probing vectors ($s \ll n_2$) to estimate the diagonal elements simultaneously. This reduces the computational cost to scale linearly with the number of genotyped animals.
4. SchurAsm-Packed: Directly assembles the matrix $S_{ped} = D + \alpha A^{-1}_{22}$ using the Schur complement identity, accumulating corrections column-by-column into a packed dense format. It then performs a direct Cholesky factorization and inversion on this assembled matrix.
5. SchurAsm-Sparse: Similar to SchurAsm-Packed but utilizes a parallel sparse LDL$^\top$ factorization on the assembled matrix.

Key Results
All five strategies reproduced the dense benchmark with high numerical accuracy, achieving Spearman rank correlations of at least 0.99.

• Accuracy: WoodColl, SchurPCG-P, and the Schur assembly variants produced near-perfect agreement with the benchmark (regression slopes $\approx$ 1.0, $R^2 = 1.0$). HutchColl showed slightly higher deviations due to its stochastic nature (mean absolute difference = 0.003, Spearman correlation = 0.999), but remained unbiased.
• Memory Efficiency:
  • DenseInv required ~18.3 GB (quadratic scaling).
  • WoodColl and HutchColl achieved the lowest memory footprints (~3.5 GB), representing an ~80% reduction. Their memory usage scales linearly with the number of genotyped animals.
  • SchurPCG-P required ~4.9 GB (73% reduction), limited by the fill-in of the sparse factorization of $A_{11}$.
  • SchurAsm-Packed required ~5.6 GB (69% reduction), scaling quadratically but with a smaller constant due to packed storage.
  • SchurAsm-Sparse was the least efficient, requiring ~30 GB (63% more than DenseInv) because sparse index overhead is detrimental when applied to a fully dense matrix.
• Computation Time:
  • DenseInv was the fastest (~3.4 min mean).
  • HutchColl was the fastest proposed method (~4.3 min mean), comparable to the dense benchmark at this scale.
  • WoodColl (66 min), SchurPCG-P (166 min), and the assembly variants (~166–194 min) were significantly slower, primarily due to the high cost of $n_2$ iterative solves or repeated factorization steps.

Significance and Claims
The paper claims that these strategies successfully extend the computation of approximate reliabilities in ssGBLUP beyond the limits of dense-inverse approaches. The authors posit that the choice of method depends on the population size and available resources:

• For populations up to ~50,000 genotyped animals, DenseInv remains feasible and fastest if memory permits.
• For intermediate populations (50,000–200,000), SchurAsm-Packed is recommended for exact results without separate $A^{-1}_{22}$ allocation, while WoodColl offers exact results with linear memory scaling.
• For large populations (>200,000), HutchColl is identified as the only practical option evaluated, as it maintains linear scaling in both memory and computation time, making it viable for populations in the hundreds of thousands or millions.

The study concludes that while WoodColl and HutchColl offer the most robust scalability by avoiding the $A_{11}$ factorization bottleneck entirely, HutchColl specifically provides a unique balance of speed and memory efficiency suitable for the largest-scale genomic evaluations. The authors note that these methods are orthogonal to APY core-size reductions for the genomic component, suggesting they can be combined to address the two main bottlenecks in ssGBLUP reliability computation.