# ASReml-R: Storing A inverse as a sparse matrix

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**Gregor Gorjanc (gg)**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I was testing ASReml-R program (an R package that links propriety ASReml binaries that can be used only with valid licence) this week and had to do some manipulations with the numerator relationship matrix (A). ASReml-R provides a function (asreml.Ainverse) that can create inverse of A directly from the pedigree as this inverse is needed in pedigree based mixed model. Bulding inverse of A directly from a pedigree is a well known result dating back to Henderson in 1970s or so. The funny thing is that it is cheaper to setup inverse of A directly than to setup up first A and then to invert it. In addition, inverse of A is very spare so it is easy/cheap to store it. Documentation for asreml.Ainverse has a tiny example of usage. Since the result of this function is a list with several elements (data.frame with “triplets” for non-zero elements of inverse of A, inbreeding coefficients, …) example also shows how to create a matrix object in R as shown bellow:

library(package="asreml") ## Create test pedigree ped

So the inverse of A would be:

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 5 0 0 -2 -2 0 0.0 0.0 0.000000 0.000000 [2,] 0 3 0 0 0 -2 0.0 0.0 0.000000 0.000000 [3,] 0 0 1 0 0 0 0.0 0.0 0.000000 0.000000 [4,] -2 0 0 3 0 1 -2.0 0.0 0.000000 0.000000 [5,] -2 0 0 0 3 1 0.0 -2.0 0.000000 0.000000 [6,] 0 -2 0 1 1 4 -2.0 -2.0 0.000000 0.000000 [7,] 0 0 0 -2 0 -2 4.5 0.5 -1.000000 0.000000 [8,] 0 0 0 0 -2 -2 0.5 4.5 -1.000000 0.000000 [9,] 0 0 0 0 0 0 -1.0 -1.0 4.909091 -2.909091 [10,] 0 0 0 0 0 0 0.0 0.0 -2.909091 2.909091

However, this is problematic as it creates a dense matrix – zero values are also stored (you can see them). If we would have 1,000 individuals, such a matrix would consume 7.6 Mb of RAM (= (((1000 * (1000 + 1)) / 2) * 16) / 2^20). This is not a lot, but with 10,000 individuals we would already need 763 Mb of RAM, which can create some bottlenecks. A solution is to create a sparse matrix using the Matrix R package. Luckily we have all the ingredients prepared by asreml.Ainverse function – the triplets. However, the essential R code is a bit daunting and I had to test several options before I figured it out – code from my previous post helped;)

## Load package library(package="Matrix") ## Number of pedigree members nI

And the inverse of A is now:

10 x 10 sparse Matrix of class "dsCMatrix" [[ suppressing 10 column names ‘1’, ‘2’, ‘3’ ... ]] 1 5 . . -2 -2 . . . . . 2 . 3 . . . -2 . . . . 3 . . 1 . . . . . . . 4 -2 . . 3 . 1 -2.0 . . . 5 -2 . . . 3 1 . -2.0 . . 6 . -2 . 1 1 4 -2.0 -2.0 . . 7 . . . -2 . -2 4.5 0.5 -1.000000 . 8 . . . . -2 -2 0.5 4.5 -1.000000 . 9 . . . . . . -1.0 -1.0 4.909091 -2.909091 10 . . . . . . . . -2.909091 2.909091

you can clearly see the structure and it soon becomes obvious why such a storage is more efficient.

If we want to go back from matrix to triplet form (this might be useful if we want to create a matrix for programs as ASReml) we can use the following code:

## Convert back to triplet form - first the matrix tmp2

To

**leave a comment**for the author, please follow the link and comment on their blog:**Gregor Gorjanc (gg)**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.