Colleau (1992) realized that computations associated with matrix the product Ax or similar products (in different notations) TDT'x , T'x , Tx etc, where T is a triangular matrix that in the i-th row contains the fractions of individual i that come from all preceding individuals, could be easily done by solving instead of computing the whole matrix A. As done for instance in Aguilar et al. 2011 ; this is in fact a key algorithm in implementation of ssGBLUP.
Whereas matrix T is a bit complicated to compute, it was realized (in the 70's ?) that $latex T^{-1} $latex had a very simple structure: $latex T^{-1} = I -P $latex where P contains 0.5 in the (individual, sire) and (individual, dam) locations and 0 otherwise - see Quaas (1988) for explanations.
A word on notation: Colleau calls T the matrix that links animals to parents; Quaas calls it P. In my opinion Quaas notation is more popular and I stick to it.
The matrix of relationship can be shown to be
$latex A = (I - P)^{- 1}D{(I - P)^{- 1}}^{'} $latex
With inverse
$latex A^{- 1} = (I - P)^{'}D^{- 1}(I - P) $latex
Where $latex P $latex contains, in the i-th row, values of 0.5 for sire(i) and dam(i). This represents the expression $latex u_{i} = u_{s} + u_{d} + \phi $latex. Matrix $latex D $latex is diagonal and contains the variance of the mendelian sampling, i.e. $latex D_{ii} = 0.5 -0.25 \left( F_{s} + F_{d} \right) $latex with Meuwissen and Luo (1992) presenting, in the case of unknown ancestor, $latex s = 0 $latex (or $latex d = 0 $latex ), their programming set $latex F_{0} = - 1 $latex.
The idea of Colleau is to use this decomposition to compute quickly products of the form $latex x = Av $latex solving the system of equations $latex A^{- 1}x = v $latex as $latex (I - P)^{'}D^{- 1}(I - P)x = v $latex from left to right in three steps:
Also, for each row, the locations of the $latex - 0.5 $latex are simply its sire and dam. For instance, the solution of $latex a_{2} $latex is $latex a_{2} = v_{2} + 0.5(a_{3} + a_{4} + a_{6}) $latex, in other words, the value of $latex v_{2} $latex plus 0.5 times the values of its offspring. Then, we can solve from the bottom to the top and adding "contributions" of +0.5 the value of "a" to the ancestors of animal i , and once we are in row i there is no more modifications to make to animal i **if ancestors are coded before offspring**. As follows:
do i=nanim,1,-1
a(i)=a(i)+v(i)
if (s(i)>0) a(s(i))=a(s(i))+0.5*a(i)
if (d(i)>0) a(d(i))=a(d(i))+0.5*a(i)
enddo
To solve (2) there are two options, one is to precompute the diagonal of D and then
b=0
b=D*a
do i=1,nanim
dii=0.5-0.25*(F(s(i)+F(d(i))
b(i)=a(i)*dii
enddo
$latex I - P = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ - 0.5 & - 0.5 & 1 & 0 & 0 & 0 \\ - 0.5 & - 0.5 & 0 & 1 & 0 & 0 \\ - 0.5 & 0 & - 0.5 & 0 & 1 & 0 \\ 0 & - 0.5 & 0 & 0 & - 0.5 & 1 \\ \end{pmatrix} $latex
To solve the system
Such that we can solve from the "top" to the "bottom".
x=0
do i=1,nanim
x(i)=b(i)
if(s(i)/=0) then
x(i)=b(i)+0.5*x(s(i))
endif
if(d(i)/=0) then
x(i)=b(i)+0.5*x(d(i))
endif
enddo
