Invited review: efficient computation strategies in genomic selection

Abstract

The purpose of this study is review and evaluation of computing methods used in genomic selection for animal breeding. Commonly used models include SNP BLUP with extensions (BayesA, etc), genomic BLUP (GBLUP) and single-step GBLUP (ssGBLUP). These models are applied for genomewide association studies (GWAS), genomic prediction and parameter estimation. Solving methods include finite Cholesky decomposition possibly with a sparse implementation, and iterative Gauss–Seidel (GS) or preconditioned conjugate gradient (PCG), the last two methods possibly with iteration on data. Details are provided that can drastically decrease some computations. For SNP BLUP especially with sampling and large number of SNP, the only choice is GS with iteration on data and adjustment of residuals. If only solutions are required, PCG by iteration on data is a clear choice. A genomic relationship matrix (GRM) has limited dimensionality due to small effective population size, resulting in infinite number of generalized inverses of GRM for large genotyped populations. A specific inverse called APY requires only a small fraction of GRM, is sparse and can be computed and stored at a low cost for millions of animals. With APY inverse and PCG iteration, GBLUP and ssGBLUP can be applied to any population. Both tools can be applied to GWAS. When the system of equations is sparse but contains dense blocks, a recently developed package for sparse Cholesky decomposition and sparse inversion called YAMS has greatly improved performance over packages where such blocks were treated as sparse. With YAMS, GREML and possibly single-step GREML can be applied to populations with >50 000 genotyped animals. From a computational perspective, genomic selection is becoming a mature methodology.

Implications

Tools used in genomic selection range from those based on estimation of SNP effects to animal models using a genomic relationship matrix (GRM). When many animals are ungenotyped, modeling options include multi-step or single-step methods. Computations involving mixed models with SNP effects or the inverse of GRM have good computing properties, whereas those including their combination or plain GRM do not. Recent discovery of efficient inverse of GRM removed computing limits from genomic analyses involving GRM. Same discovery enables GREML for very large data sets.

Introduction

The concept of genomic selection (Meuwissen et al., 2001) generated great excitement in the animal breeding community. With genomic information from SNP panels, one can achieve accuracy from young animals almost as high as from a progeny test (Schaeffer, 2006; VanRaden et al., 2009) at a greatly reduced cost. Initially the genomic computations used versions of BLUP with SNP fitted as random effects, such as SNP BLUP, BayesB, etc., (see Gianola et al., 2009 for a review). An alternative form of SNP BLUP is genomic BLUP (GBLUP), where the animal effect is fit with a GRM (VanRaden, 2008). For prediction, SNP BLUP and GBLUP are equivalent models (VanRaden, 2008; Strandén and Christensen, 2011).

When only a small fraction of the population is genotyped, the information from ungenotyped animals can be summarized in pseudo-observations for genotyped animals. Alternatively, GBLUP can be extended to single-step GBLUP (ssGBLUP) (Aguilar et al., 2010; Christensen and Lund, 2010), where a numerator relationship matrix (NRM) for all individuals and GRM are combined and then applied to BLUP. Benefits of ssGBLUP include simplicity of application (another BLUP), avoidance of double counting, and accounting for pre-selection on Mendelian sampling (Legarra et al., 2014). ssGBLUP was extended to be based partially or completely on SNP effects (Fernando et al., 2014; Liu et al., 2014).

Large interest in the genomic community exists in identifying an optimal set of SNP and their variances for increased accuracy of evaluation. Such an identification is straightforward with Bayesian SNP models (e.g. Gianola et al., 2009). As SNP BLUP and GBLUP are equivalent, SNP and their variances selected in Bayesian SNP models can be transformed to weighted GRM and subsequently BLUP. Methods exist to generate SNP weights directly in the GBLUP model (VanRaden, 2008; Zhang et al., 2010; Sun et al., 2012) or ssGBLUP (Wang et al., 2012).

Initially the extent of genomic selection was limited by costs of genotyping. With these costs dropping dramatically, the number of genotyped animals has increased dramatically, with very high-density SNP chips available. In dairy cattle, over 1 million Holsteins had been genotyped in the United States as of March 2016 (Council on Dairy Cattle Breeding; https://www.cdcb.us/Genotype/cur_density.html). Subsequently, computing costs are becoming more of an issue. The purpose of this paper is to present and discuss computational methods in genomic selection.

Definition of terms and basic models

Let A denote a NRM based on pedigrees. Let Z be a matrix of gene content with z _ij being the number of occurrences of the minor allele in SNP j of animal i. The values of ‘raw’ z _ij are 0, 1 and 2, whereas that of centered z _ij is −2p _j , 1−p _j and 2−p _j , where p _j is gene frequency of SNP j. The simplest SNP BLUP model is

$$y_{i} \,{\equals}\,\mu {\plus}\mathop{\sum}\limits_j^{} {z_{{ij}} a_{j} {\plus}e_{i} } $$

where y _i is a (pseudo)phenotype of animal i, $${\bf a}\,\sim\,N\left( {0,{\bf I}\sigma _{a}^{2} } \right)$$ the vector of SNP effects, $${\bf e}\,\sim\,N\left( {0,{\bf I}\sigma _{e}^{2} } \right)$$ the vector of residuals, and $$\sigma _{a}^{2} $$ and $$\sigma _{e}^{2} $$ are SNP and residual variances, respectively. GBLUP can be defined as

$$y_{i} \,{\equals}\,\mu \,{\plus}\,u_{i} \,{\plus}\,e_{i} $$

where $${\bf u}\,\sim\,N\left( {0,{\bf G}\sigma _{u}^{2} } \right)$$ is a vector of animal effects, G the GRM and $$\sigma _{u}^{2} $$ the additive variance. G can be derived by a transformation from SNP BLUP:

$${\bf{G}} \,{\equals}\,{\bf{ZZ'}} {{\sigma _{a}^{{\rm 2}} } \over {\sigma _{u}^{{\rm 2}} }}$$

possibly with an alternative scaling factor (VanRaden, 2008):

$${\bf{G}} \,{\equals}\,{{{\bf{ZZ'}} } \over {\mathop{\sum}\limits_j^{} {\left[ {p_{j} \left( {{\rm 1}{\minus}p_{j} } \right)} \right]} }}$$

If SNP BLUP involves heterogeneous variance for SNP: $${\bf a}\,\sim\,N(0,{\bf D}\sigma _{a}^{2} )$$ , where D is a diagonal matrix of weights, the equivalent ‘weighted’ G is

$${\bf{G}} \,{\equals}\,{{{\bf{ZDZ'}} } \over {\mathop{\sum}\limits_j^{} {\left[ {p_{j} \left( {{\rm 1}{\minus}p_{j} } \right)} \right]} }}$$

SNP BLUP and GBLUP are equivalent models (Strandén and Christensen, 2011).

When only a fraction of animals is genotyped, the numerator and genomic relationships can be combined (Legarra et al., 2009; Christensen and Lund, 2010):

$${\bf H}\,{\equals}\;{\bf A}{\plus}\left[ {\matrix{ {{\bf A}_{{12}} {\bf A}_{{22}}^{{{\minus}1}} } & {\bf 0} \cr {\bf 0} & {\bf I} \cr } } \right]\left[ {\matrix{ {\bf I} \cr {\bf I} \cr } } \right]({\bf G}{\minus}{\bf A}_{{22}} )\left[ {\matrix{ {\bf I} & {\bf I} \cr } } \right]\left[ {\matrix{ {{\bf A}_{{22}}^{{{\minus}1}} {\bf A}_{{21}} } & {\bf 0} \cr {\bf 0} & {\bf I} \cr } } \right]$$

where H is a combined matrix and indices in A refer to ungenotyped (1) and genotyped (2) animals. The inverse of matrix H is (Aguilar et al., 2010; Christensen and Lund, 2010)

$${\bf H}^{{{\minus}1}} \,{\equals}\;{\bf A}^{{{\minus}1}} {\plus}\left[ {\matrix{ {\bf 0} & {\bf 0} \cr {\bf 0} & {{\bf G}^{{{\minus}1}} {\minus}{\bf A}_{{22}}^{{{\minus}1}} } \cr } } \right]$$

BLUP with matrix H is called ssGBLUP.

Computing operations in animals breeding

Common operations in animal breeding include solving of mixed model equations and variance component estimation either via REML or Bayesian methods with Gibbs sampling. Such operations in the animal breeding context are described in Mrode (2014) and Misztal (2014), and are typically based on Henderson’s mixed model equations (MME). Let this linear system of equations be

$${\bf{B}} \,{\bf{x}} \,{\bf{{\equals}}} \;{\bf{y}} $$

where B is left-hand-side (LHS), y the right-hand side, and x the vector of solutions. Solving mixed model equations can be accomplished by finite and iterative methods. In finite methods, exact solutions (except for numerical errors) are obtained in a finite number of steps. Popular methods for general matrices are based on the LU decomposition (B=LU, where L and U are lower and upper triangular, respectively) and for symmetric matrices the Cholesky decomposition (B=LL', where L is lower triangular). As often mixed model equations are not full rank, finite methods used for mixed models need to compute generalized solutions rather than cause numerical exceptions. For example, the Cholesky decomposition can be modified to

$${\bf{B\,{\equals}\,LDL}} \prime$$

where D is a diagonal matrix with zero elements for redundant equations and non-zeros otherwise. Although finite methods usually compute solutions with high precision, their cost is high, usually close to cubic. Finite methods can be used to calculate inverses, which are used in REML or for prediction error variance (PEV) calculations. In general, the cost of the inverse is a few times that of a solution by decomposition.

When matrices are sparse, with nearly all elements of LHS equal to 0, a ‘sparse’ Cholesky decomposition exists that avoids operations with zeros, leading to about quadratic rather than cubic costs (George and Liu, 1981). Also, an algorithm exists for computing a ‘sparse’ inverse, where elements computed are those corresponding to non-zeroes in the original matrix (Takahashi et al., 1973). Such an inverse is sufficient to calculate traces in REML (Misztal and Perez-Enciso, 1993) or PEV. However, equations involving SNP effects are usually dense. When the system of equations is sparse but contains dense blocks, a sparse matrix package can be modified for greatly improved performance in REML and other applications (Masuda et al., 2015). This is quite new as most matrix packages are efficient for either dense or sparse matrices, but not for mixtures of both.

Iterative methods are such that progressively more accurate solutions are obtained every next round. Two such methods are popular in animal breeding. The first one is Gauss–Seidel iteration (GS), which is also a backbone of Gibbs sampling. GS is very simple, but requires access to equations by rows. The second one is preconditioned conjugate gradient (PCG), a method difficult to understand, but easy to implement and usually with superior convergence rate to GS (Tsuruta et al., 2001).

In animal breeding application, the size of LHS can be much larger than the size of data to generate LHS. For example, one line of data in a 10 trait model with 10 effects per trait includes 20 numbers (one for each trait and effect) but generates 10 000 contributions to LHS. A special implementation of iteration algorithms is matrix-free or by iteration on data (Schaeffer and Kennedy, 1986; Misztal and Gianola, 1987), where coefficients of LHS are regenerated from the data every round of iteration. Although the implementation of iteration on data by GS can be complicated as the coefficients need to be recreated row by row, such an implementation in PCG is trivial as PCG only requires a product of B by a vector (Bq, with q a given vector), with no individual elements necessary (except for diagonals).

Tricks and rules in mixed model and genomic computations

Product of numerator relationship matrix

Computing the NRM A by a tabular method has a quadratic cost. Also, computing relationships for specific animals requires creating a matrix for all ancestors. For large populations, A would not fit into memory. When only a product A by a vector, say q, is needed, Aq can be computed inexpensively by two scans of the pedigree (Colleau, 2002). By selecting zeros in q, one can also compute a product of a section of A by a vector (e.g. A ₂₂ q) (Misztal et al., 2009; Aguilar et al., 2011).

Sequential multiplication and preconditioned conjugate gradient

The LHS of the mixed model equations contains expressions like Z'Z. The PCG algorithm by iteration on data requires computing Z'Zq, where q is a given vector. Strandén and Lidauer (1999) found that while the cost of (Z'Z)q is very high, the cost of Z'(Zq) is much lower.

Convergence properties with various systems of equations

When mixed models are solved by iteration (GS or PCG), the convergence rate is better with fewer equations per phenotype. Experiences indicate that the convergence rate is very good with sire models, much slower with animal models, and could be especially slow with young animals not tied to phenotypes or progeny.

The regular MME equations involve an inverse of A or G, e.g.

$$\left( {{\bf{I}} {\plus}{\bf{A}} ^{{{\minus}1}} \alpha } \right) {\bf{ {\vskip -3pt \hskip 2pt {\hat }}{\hskip -2.5pt a}\,{\equals}\;y}} $$

An equivalent system of equation avoids creating an inverse (Strandén and Garrick, 2009):

$$\left( {{\bf{A{\plus}I}} \alpha } \right){\bf{\bf{ {\vskip -3pt \hskip 2pt {\hat }}{\hskip -2.5pt a}\;{\equals}\;Ay}} $$

The convergence rate in such models is poor and may not occur for larger systems of equations. This is because A ⁻¹ ‘connects’ only parents and progeny resulting in ‘sparse’ equations, whereas all animals in A are connected, resulting in dense equations. As G differs little from A (VanRaden, 2008; Wang et al., 2014), the same thinking applies to models with G and G ⁻¹.

Computations in SNP BLUP (and BayesX)

SNP BLUP and its various forms (including BayesA, BayesB, etc.) are very popular models in genomic computations. When the number of SNP is small, the MME for SNP BLUP can be created explicitly. As the memory and compute requirements are quadratic, explicit storage of SNP BLUP equations is unfeasible with large number of SNP. For example, memory requirements for MME with 1 million SNP would be 10¹² elements also with 10¹² elements of MME created from each genotyped animal. Various options in such a case have been explored by Legarra and Misztal (2008) with a data set on about 2 k individuals and 20 k SNP markers. Solution by Cholesky decomposition took 2 h, and computations would increase cubically with the number of SNP, whereas memory would increase quadratically. With solution by GS and straightforward iteration on data, computations took about 4 days, with linear memory and quadratic memory costs. With a special GS implementation called GSRU, where currently computed solutions are used to adjust residuals for all individuals, computing took about 1 min and both computing and memory costs were linear. Finally, with PCG iteration on data, computing took only 20 s, with approximately linear costs for memory and computations. Although GSRU is the only realistic algorithm with millions of SNP for Bayesian models, the PCG approach can be useful for models where SNP variances are derived from SNP solutions (Zhang et al., 2010; Sun et al., 2012; Wang et al., 2012). A hybrid method would use PCG iteration with inner (non-Monte Carlo-based) updates of weights (VanRaden, 2008).

Variations of single-step equations

Due to the prohibitive cubic cost of inversion of G with many genotyped animals, alternative equations were proposed for ssGBLUP. Equations in Misztal et al. (2009) used H rather than H ⁻¹. Convergence was achieved only with a small number of animals. Equations by Legarra and Ducrocq (2012) used G and A ₂₂, which are easy to use, instead of G ⁻¹ and A ₂₂ ⁻¹. The convergence rate was slow with medium models and not achieved with large models (Aguilar et al., 2013). Liu et al. (2014) proposed equations with SNP effects for genotyped animals. They reported lack of convergence with real data. All these models illustrate computing advantages with inverses of dispersion matrices (A ⁻¹ or G ⁻¹). As SNP BLUP models have good computing properties, Fernando et al. (2014) looked into ssGBLUP where ungenotyped animals were imputed and only SNP effects (plus imputation errors for ungenotyped animals) were estimated. Although convergence rates were good, imputations required a supercomputer and imputed genotypes required large memory for big populations. Better computing times were reported using a computer with fast graphic cards (Golden et al., 2016). Meuwissen et al. (2015) presented a new type of H matrix based on (possibly fractional) imputation of ancestors by segregation analyses, followed by a sophisticated tailoring of G to A to compensate for errors in imputation. Although the cost of segregation analysis is high, it is not yet clear whether such an approach offers accuracy benefits when most ancestors are genotyped, and the process of tailoring is rather complex. All non-traditional ssGBLUP would require new research and investment in code for implementation of more complex models, such as multiple-trait, random regression or maternal models.

Sparse genomic relationship matrix

Computations with ssGBLUP would be easy for large number of animals if G ⁻¹ was sparse, like A ⁻¹. Past studies suggested that G had a limited rank, as inversion is unstable with >5 to 10 k animals, and Gs as used in genetic evaluations are blended with A ₂₂ (VanRaden, 2008; Aguilar et al. 2010) to achieve full rank. The eigenvalue decomposition of G is:

$${\bf{G}} \,\,{\bf{{\equals}}} \;{\bf{UDU}} '$$

where U is a matrix of eigenvectors and D the matrix of eigenvalues. If all eigenvalues are positive, the inverse of G is

$${\bf{G}} ^{{{\minus}1}} \,\,{\bf{{\equals}}} \;{\bf{U}} '{\bf{D}} ^{{{\minus}1}} {\bf{U}} $$

If some eigenvalues are 0, the inverse does not exist. If, due to noise, many eigenvalues are close to 0 and carry no information, they can cause instability of the inverse through large elements of D ⁻¹.

Let D _t indicate a fraction of D with non-negligible eigenvalues, and let U _t be corresponding eigenvalues. Then

$${\bf{G}} ^{{\minus}} \,{\bf{{\equals}\;U}} _{t}^{'} {\bf{D}} _{t}^{{{\rm {\minus}1}}} {\bf{U}} _{t} $$

If G ⁻¹ is to be created explicitly, the computing cost is cubic with quadratic storage. If G ⁻ is used by a PCG algorithm as

$${\bf{G}} ^{{\minus}} \,{\bf{q\;{\equals}}}\, \left( {{\bf{U}} _{t}^{'} \left( {{\bf{D}} _{t}^{{{\rm {\minus}1}}} \left( {{\bf U}_{t}^{{}} {\bf q}} \right)} \right)} \right)$$

and the number of eigenvalues in D _t is small, G ⁻ q can be calculated and stored at linear cost with respect to the number of genotyped animals. However, computing U and D is expensive and requires storage of G in full.

When G has a limited dimensionality, Misztal et al. (2014) proposed to calculate G ⁻¹ using only a fraction of G and without eigenvalue decomposition. Let

$${\bf{G\,{\equals}}}\, \left[ {\matrix{ {{\bf{G}} _{{cc}} } & {{\bf{G}} _{{cn}} } \cr {{\bf{G}} _{{cn}} } & {{\bf{G}} _{{nn}} } \cr } } \right]$$

where index c denotes core animals and index n non-core animals. In the algorithm for proven and young animals (APY)

$${\bf G}^{{{\minus}1}} \,\approx\,\left[ {\matrix{ {{\bf G}_{{cc}}^{{{\minus}1}} } & {\bf 0} \cr {\bf 0} & {\bf 0} \cr } } \right]{\plus}\left[ {\matrix{ {{\minus}{\bf G}_{{cc}}^{{{\minus}1}} {\bf G}_{{cn}} } \cr {\bf I} \cr } } \right]{\bf M}^{{{\minus}1}} \left[ {\matrix{ {{\minus}{\bf G}_{{nc}} {\bf G}_{{cc}}^{{{\minus}1}} } & {\bf I} \cr } } \right],$$

where M is a diagonal matrix with elements $m_{i} \,{\equals}\,g_{{ii}} {\minus}{\bf G}_{{ic}} {\bf G}_{{cc}}^{{{\minus}1}} {\bf G}_{{ci}} $ for individual i in the core group. Inversion of G by APY has a cubic cost (and quadratic memory) for core individuals and a linear cost (and linear memory) for non-core individuals. Figure 1 shows an example non-zero structure of a regular inverse (after blending to make it positive-definite) and different APY inverses. The inverses are diagonal for non-core animals. As the regular G is singular for a large population, an infinite number of generalized inverses exist, including those by the APY algorithm with different sets of core animals.

Figure 1 Non-zero pattern of a regular (a) and APY (b to d) inverses of the genomic relationship matrix with different choices of core animals.

Core animals can be randomly selected (Fragomeni et al., 2015), with an optimal number of about 12 000 for Holsteins (Masuda et al., 2016) and 10 000 for Angus (Lourenco et al., 2015). The number of core animals is linked to effective population size (Misztal, 2016). Optimal number of core animals is equal to the number of the largest eigenvalues in D explaining 98% of the variance in G (Pocrnic et al., 2016a). With such a number, prediction of breeding values with APY was more accurate than with a regular inverse. Computation of G ⁻¹ for 570 k genotyped Holsteins by APY took only about 2 h, whereas a month computing with a large memory would be required for a regular inverse (Masuda et al., 2016). The APY algorithm removes computing limits from inversion of G and subsequently ssGBLUP.

Computing of A₂₂

ssGBLUP as in Aguilar et al. (2010) requires $${\bf A}_{{22}}^{{{\minus}1}} $$ . Although A ⁻¹ is sparse, $${\bf A}_{{22}}^{{{\minus}1}} $$ may be relatively dense (Faux and Gengler, 2013), and its storage for potentially millions of genotyped animals impossible. This matrix can be calculated indirectly (Strandén and Mäntysaari, 2014).

$${\bf{A}} _{{22}}^{{{\minus}1}} \,{\bf{{\equals}\,A}} ^{{22}} {\minus}\left( {{\bf{A}} ^{{12}} } \right)'\left( {{\bf{A}} ^{{11}} } \right)^{{{\minus}1}} {\bf{A}} ^{{12}} $$

where all submatrices are sparse and can be stored explicitly. With PCG iteration, explicit $$A_{{22}}^{{{\minus}1}} $$ is not required, and its product with a vector q can be calculated every round as follows (Masuda et al., 2016):

$${\bf{A}} _{{22}}^{{{\minus}1}} {\bf{q\,{\equals}}}\, \left[ {{\bf{A}} ^{{22}} {\minus}\left( {{\bf{A}} ^{{12}} } \right)'\left( {{\bf{A}} ^{{11}} } \right)^{{{\minus}1}} {\bf{A}} ^{{12}} } \right]{\bf{q}} $$

In particular

$${\mib{s}} \,{\bf{{\equals}\, }} \left( {{\bf{A}} ^{{11}} } \right)^{{{\minus}1}} {\bf{A}} ^{{12}} {\bf{q}} $$

is computed as a solution to:

$${\bf{A}} ^{{11}} {\mib{s}} \,\,{\bf{{\equals}\,\, A}} ^{{12}} {\bf{q}} $$

As all matrices are sparse, the cost of above calculations is small (Masuda et al., 2016).

GREML and APY

GREML is REML using a GRM. Similarly, REML using matrix H can be called single-step GREML (ssGREML). In the past, the success with REML in animal populations was due to sparse inverse of A ⁻¹ and sparse matrix factorization/inversion (Misztal and Perez-Enciso, 1993; Pérez-Enciso et al., 1994). If G is dense and large, Masuda et al. (2015) showed that GREML with the old sparse matrix package is inefficient and unstable. A new package called YAMS recognizes dense blocks in MME, rearranges computations accordingly, allowing also for parallel computations for large dense blocks. With YAMS, GREML was much faster and more reliable.

Additional issues

This paper focused on algorithmic issues relevant to current mixed model computations in genomics in animal breeding. Many other issues exist and many will become relevant. For example, genotyping of animals with small amount of information leads to greater importance of avoidance of double counting in multi-step evaluations, or of more accurate match between genomic and additive relationships in ssGBLUP. With efforts to use sequence data, important operations are aligning short fragments, imputation, genomewide association studies, and incorporations of the detected variants in the evaluation. In dairy, a de facto global breeding scheme for Holsteins creates a need for a comprehensive genomic evaluation by Interbull that has a high accuracy, yet does not compromise privacy of the member countries. In general, while brute-force approaches may be formidable, as or more accurate yet inexpensive approaches may exist that exploit peculiarities of data including a limited effective population in farm populations. For instance, the sequence includes three G base pairs, but farm populations can be described by about 4000 to 20 000 chromosomal segments (Pocrnic et al., 2016b). This means a block size of about 150 to 750 kbases with hard to differentiate SNPs, and ability to reduce costs by working with blocks rather than individual SNP. This also suggests that increases of genomic accuracies past 4000 to 20 000 of high accuracy genotyped animals are limited (ignoring G×E and decays of predictivity over time). Similarly, if Interbull countries are unwilling to share genotypes on individual bulls or SNP effects derived from their population, they can contribute GEBV and genotypes of artificially generated animals (Maantysaari, 2015, personal communication) with their number and accuracies relevant to each country.

Conclusions

In conclusion, common operations in genetic selection are computationally feasible for an almost unlimited number of genotyped individuals using appropriate algorithms. Solving ‘SNP’ equations is best accomplished by iteration on data using either a special form of GS iteration or the PCG algorithm. Solving ssGBLUP is best accomplished by iteration on data using the PCG algorithm. For GREML, a sparse matrix package needs to account for dense blocks. Computations in the two latter cases greatly benefit from efficient (and sparse) inverse of the GRM.