Computing A and its inverse

Read pedigree data
Computing A matrix
Computing the inverse of A matrix

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
 )