make_invA

makeA(pedigree)

Constructs the inverse of Wright's Numerator Relationship Matrix (WNRM) directly from a given pedigree structure. This will be more efficient that computing the inverse in NumPy using inv(sigma) because it exploits structural properties.

Parameters

• pedigree : dict
  Pedigree structure in an {int: [int, int]} dictionary format, such as that returned from load_ped.

Returns

• ndarray
  A dense NumPy array which is equal to the inverse of Wright's Numerator Relationship Matrix.

Examples

Consider an example cohort with four candidates, which when loaded into Python produces the following dictionary representing the pedigree structure

>>> pedigree = robustocs.load_ped("example.ped")
>>> pedigree
{1: [0, 0], 2: [1, 0], 3: [1, 2], 4: [1, 3]}

We can then use this in make_invA to form the WNRM ($A$) as

>>> robustocs.make_invA(pedigree)
array([[ 2.4047619 , -0.16666667, -0.42857143, -1.14285714],
       [-0.16666667,  1.83333333, -1.        ,  0.        ],
       [-0.42857143, -1.        ,  2.57142857, -1.14285714],
       [-1.14285714,  0.        , -1.14285714,  2.28571429]])

which we see is identical to the matrix produced by

>>> numpy.linalg.inv(A)
array([[ 2.4047619 , -0.16666667, -0.42857143, -1.14285714],
       [-0.16666667,  1.83333333, -1.        ,  0.        ],
       [-0.42857143, -1.        ,  2.57142857, -1.14285714],
       [-1.14285714,  0.        , -1.14285714,  2.28571429]])

only computed faster due to exploiting the pedigree structure [1].

To see an example of this discrepancy in speed, we look at a slightly larger pedigree

i,p,q
1,0,0
2,0,0
3,1,0
4,1,2
5,3,4
6,1,4
7,5,6

and compute the inverse 10000 times with each method. We find that it takes make_invA(ped) 0.59 seconds and inv(makeA(ped)) 0.99 seconds, with both forming the inverse

 [[ 2.333  0.5   -0.667 -0.5    0.    -1.     0.   ]
 [ 0.5    1.5    0.    -1.     0.     0.     0.   ]
 [-0.667  0.     1.833  0.5   -1.     0.     0.   ]
 [-0.5   -1.     0.5    3.    -1.    -1.     0.   ]
 [ 0.     0.    -1.    -1.     2.615  0.615 -1.231]
 [-1.     0.     0.    -1.     0.615  2.615 -1.231]
 [ 0.     0.     0.     0.    -1.231 -1.231  2.462]]

References

• [1] Henderson (1976), A Simple Method for Computing the Inverse of a Numerator Relationship Matrix Used in Prediction of Breeding Values, Biometrics 32:1, pg. 69-83.