Computing A and its inverse

We will evaluate what is the fastest way to compute an \(\mathbf{A}\) matrix and its inverse.

Read pedigree data

library(readxl)
dat <- read_excel("../data/Simdata.xlsx")

Computing A matrix

Tabular method

The createA function computes the \(\mathbf{A}\) matrix using a tabular method described in Wright, 1922, Cruden, 1949, and Emik and Terrill, 1949.

args(createA)
function (s, d) 
NULL
A_tab <- createA(s=dat$Sire, d=dat$Dam)

Recursive method for computing L

The createL function computes the \(\mathbf{L}\) matrix described in Henderson, 1976. Then, \(\mathbf{A}\) can be obtained as \(\mathbf{A} = \mathbf{LL}'\).

args(createL)
function (s, d) 
NULL
L <- createL(s=dat$Sire, d=dat$Dam)
A_L <- L%*%t(L)

Are the two matrices equal?

all.equal(A_tab, A_L)

Which algorithm is the fastest to compute A?

The microbenchmark function in the microbenchmark package can be used to benchmark the efficiency of your R code. We will use this function to compare the performance of tabular method (createA) and recursive method for computing L (createL).

library(microbenchmark) 
microbenchmark(
   "A = tabular method" = {
     A_tab <- createA(s=dat$Sire, d=dat$Dam)
     },
   
   "A = LL'" = {
     L <- createL(s=dat$Sire, d=dat$Dam)
     A_L <- L%*%t(L) 
     },
   
   times = 50L
 ) 

Alternative appraoch is to replace L%*%t(L) with tcrossprod(L).

microbenchmark(
   "A = tabular method" = {
     A_tab <- createA(s=dat$Sire, d=dat$Dam)
     },
   
   "A = LL'" = {
     L <- createL(s=dat$Sire, d=dat$Dam)
     A_L <- tcrossprod(L) 
     },
   
   times = 50L
 ) 

Computing the inverse of A matrix

A simple method using L

This was first proposed in Henderson, 1976.

args(createAinvL)
function (s, d) 
NULL
Ainv1 <- createAinvL(s=dat$Sire, d=dat$Dam)

Rapid computation of the diagonal of L

This was first proposed in Quaas, 1976.

args(quaas)
function (s, d) 
NULL
Ainv2 <- quaas(s=dat$Sire, d=dat$Dam)

A simple method NOT using L designed for a non-inbred population

This was first proposed in Henderson, 1976.

args(createAinv)
function (s, d) 
NULL
Ainv3 <- createAinv(s=dat$Sire, d=dat$Dam)

Are the three matrices equal?

all.equal(Ainv1, Ainv2)
all.equal(Ainv1, Ainv3)
all.equal(Ainv2, Ainv3)

Which algorithm is the fastest to compute the inverse of A?

We will compare the three algorithms and the solve() function using the microbenchmark function.

library(microbenchmark) 
microbenchmark(
   "Ainv1 = Henderson (1976)" = {
     Ainv1 <- createAinvL(s=dat$Sire, d=dat$Dam)
     },
   
   "Ainv2 = Quaas (1976)" = {
     Ainv2 <- quaas(s=dat$Sire, d=dat$Dam)
     },
   
    "Ainv3 = Henderson (1976, non-inbred)" = {
     Ainv3 <- createAinv(s=dat$Sire, d=dat$Dam)
     },
   
   "Ainv4 = solve() function" = {
     Ainv4 <- solve(A_L)
     },
   times = 50L
 ) 

Gota Morota

February 1, 2018