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.
To help ideas, a small pedigree due to Kempthorne is as follows :
A 0 0
B 0 0
D A B
E A D
F B E
Z A B
Renumbered as
1 0 0
2 0 0
3 1 2
4 1 2
5 1 3
6 2 5
The corresponding matrices are
$latex P = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0.5 & 0.5 & 0 & 0 & 0 & 0 \\ 0.5 & 0.5 & 0 & 0 & 0 & 0 \\ 0.5 & 0 & 0.5 & 0 & 0 & 0 \\ 0 & 0.5 & 0 & 0 & 0.5 & 0 \\ \end{pmatrix} $latex
$latex D = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0.5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.5 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0.5 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0.4375 \\ \end{pmatrix} $latex
With $latex F = \begin{pmatrix} 0 & 0 & 0 & 0 & 0.25 & 0.125 \\ \end{pmatrix} $latex
And finally , A=
\[,1\] \[,2\] \[,3\] \[,4\] \[,5\] \[,6\]
\[1,\] 1.000 0.000 0.500 0.5 0.75 0.375
\[2,\] 0.000 1.000 0.500 0.5 0.25 0.625
\[3,\] 0.500 0.500 1.000 0.5 0.75 0.625
\[4,\] 0.500 0.500 0.500 1.0 0.50 0.500
\[5,\] 0.750 0.250 0.750 0.5 1.25 0.750
\[6,\] 0.375 0.625 0.625 0.5 0.75 1.125
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:
1. Solve $latex (I - P)^{'}a = v $latex
2. Solve $latex D^{- 1}b = a $latex
3. Solve $latex (I - P)x = b $latex
To solve (1) we use the special structure of
$latex (I - P)^{'} = \begin{pmatrix} 1 & 0 & - 0.5 & - 0.5 & - 0.5 & 0 \\ 0 & 1 & - 0.5 & - 0.5 & 0 & - 0.5 \\ 0 & 0 & 1 & 0 & - 0.5 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & - 0.5 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix} $latex
To solve the system
$latex \begin{pmatrix} 1 & 0 & - 0.5 & - 0.5 & - 0.5 & 0 \\ 0 & 1 & - 0.5 & - 0.5 & 0 & - 0.5 \\ 0 & 0 & 1 & 0 & - 0.5 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & - 0.5 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} a_{1} \\ a_{2} \\ \ldots \\ \\ \\ a_{6} \\ \end{pmatrix} = \begin{pmatrix} v_{1} \\ v_{2} \\ \ldots \\ \\ \\ v_{6} \\ \end{pmatrix} $latex
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:
a=0
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
Another option is to compute D on the run :
! previously set F(0) to 0
do i=1,nanim
dii=0.5-0.25*(F(s(i)+F(d(i))
b(i)=a(i)*dii
enddo
To solve (3) we use the special structure of
$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
$latex \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}\begin{pmatrix} x_{1} \\ x_{2} \\ \ldots \\ \\ \\ x_{6} \\ \end{pmatrix} = \begin{pmatrix} b_{1} \\ b_{2} \\ \ldots \\ \\ \\ b_{6} \\ \end{pmatrix} $latex
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
with these three steps we are done.
