Colleau's (2002) algorithm for dummies

 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 $T^{-1}$ had a very simple structure: $T^{-1}  = I -P$ 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

$A = (I - P)^{- 1}D{(I - P)^{- 1}}^{'}$

With inverse

$A^{- 1} = (I - P)^{'}D^{- 1}(I - P)$

Where $P$ contains, in the i-th row, values of 0.5 for sire(i) and dam(i). This represents the expression $u_{i} = u_{s} + u_{d} + \phi$. Matrix $D$ is diagonal and contains the variance of the mendelian sampling, i.e. $D_{ii} = 0.5 -0.25 \left( F_{s} + F_{d} \right)$ with Meuwissen and Luo (1992) presenting, in the case of unknown ancestor, $s = 0$  (or $d = 0$ ), their programming set $F_{0} = - 1$.

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

$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}$

$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}$

With $F = \begin{pmatrix} 0 & 0 & 0 & 0 & 0.25 & 0.125 \\ \end{pmatrix}$

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 $x = Av$ solving the system of equations $A^{- 1}x = v$ as $(I - P)^{'}D^{- 1}(I - P)x = v$ from left to right in three steps:

1.  Solve $(I - P)^{'}a = v$

2.  Solve $D^{- 1}b = a$

3.  Solve $(I - P)x = b$

To solve (1) we use the special structure of

$(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}$

To solve the system

$\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}$

Also, for each row, the locations of the $- 0.5$  are simply its sire and dam. For instance, the solution of $a_{2}$ is $a_{2} = v_{2} + 0.5(a_{3} + a_{4} + a_{6})$, in other words, the value of $v_{2}$  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

$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}$

To solve the system

$\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}$

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.