4. Algorithm to Invert Additive Relationship Matrix#

4.1. Additive Relationship Matrix: \(\mathbf{A}\)#

  • \(a_{ij}\) is two times the kinship probability (\(F_{ij}\)), which is the probability that randomly sampled homologous alleles from \(i\) and \(j\) are IBD

  • Homologous alleles are IBD if they can be traced to the same ancestral allele in a founder

  • Thus, two different alleles in founders are not IBD

  • Homologous alleles are randomly inherited one from each of the parents \(m_i\) and \(f_i\) of \(i\)

  • If \(j\) is not a descendant of \(i\), the genetic relationship of \(j\) with \(i\) is entirely through the parents of \(i\)

4.2. Recursive formula:#

\[ F_{ij} = \frac{1}{2}(F_{jm_i} + F_{jf_i}) \]

and

\[ F_{ii} = \frac{1}{2}(1 + F_{m_if_i}) \]

Similarly,

\[ a_{ij} = \frac{1}{2}(a_{jm_i} + a_{jf_i}) \]

and

\[ a_{ii} = \frac{1}{2}(2 + a_{m_if_i}) \]

4.3. Tabular Method to Compute \(\mathbf{A}\)#

  • Number individuals such that parents precede offspring

  • For founders, enter 1 on diagonal and 0 on off-diagonals

  • For non-founder \(i\) calculate row elements 1 to \(i-1\) as the average of the parental row elements

  • Set diagonal element \(i\) to \(\frac{1}{2}(2 + a_{m_if_i})\)

  • Fill columns by symmetry.

4.4. Tabular method in matrix notation#

\[\begin{split} \begin{align} \mathbf{A}_i &= \begin{bmatrix} \mathbf{A}_{i-1} & \mathbf{A}_{i-1}\mathbf{q}_i \\ \mathbf{q}'_i\mathbf{A}_{i-1} & \frac{1}{2}(2 + a_{m_if_i}) \end{bmatrix}\\ &= \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{q}'_i& 1 \end{bmatrix} \begin{bmatrix} \mathbf{A}_{i-1} & \mathbf{0} \\ \mathbf{0} & a_{ii} - \mathbf{q}'_i\mathbf{A}_{i-1}\mathbf{q}_i \end{bmatrix} \begin{bmatrix} \mathbf{I} & \mathbf{q}_i \\ \mathbf{0}& 1 \end{bmatrix} \\ &= \mathbf{LDL}' \end{align} \end{split}\]

where \(a_{ii} = \frac{1}{2}(2 + a_{m_if_i})=1 + \frac{a_{m_if_i}}{2}\), and \(\mathbf{q}_i\) has only two non-zero elements that are equal to 0.5 at positions corresponding to \(m_i\) and \(f_i\).

4.5. Numerical Example#

Pedigree:

\[\begin{split} \begin{matrix} i & f & m\\ 1 & 0 & 0\\ 2 & 0 & 0\\ 3 & 1 & 2\\ 4 & 1 & 2\\ 5 & 3 & 4 \end{matrix} \end{split}\]
A2 = [1 0
      0 1]

2×2 Matrix{Int64}:
 1  0
 0  1

q   = [0.5, 0.5]

2-element Vector{Float64}:
 0.5
 0.5

q'A2

1×2 adjoint(::Vector{Float64}) with eltype Float64:
 0.5  0.5

amf = A2[1,2]
A3 = [A2    A2*q
      q'A2  1 + 0.5amf   ]

3×3 Matrix{Float64}:
 1.0  0.0  0.5
 0.0  1.0  0.5
 0.5  0.5  1.0

q   = [0.5, 0.5, 0]
amf = A3[1,2]
A4 = [A3    A3*q
      q'A3  1 + 0.5amf   ]

4×4 Matrix{Float64}:
 1.0  0.0  0.5  0.5
 0.0  1.0  0.5  0.5
 0.5  0.5  1.0  0.5
 0.5  0.5  0.5  1.0

q   = [0, 0, 0.5, 0.5]
amf = A4[3,4]
A5 = [A4    A4*q
      q'A4  1 + 0.5amf   ]
A = A5

5×5 Matrix{Float64}:
 1.0  0.0  0.5   0.5   0.5
 0.0  1.0  0.5   0.5   0.5
 0.5  0.5  1.0   0.5   0.75
 0.5  0.5  0.5   1.0   0.75
 0.5  0.5  0.75  0.75  1.25

4.6. Inverse of \(\mathbf{A}_i\)#

We have seen that \(\mathbf{A}_i\) can be written as: $\( \mathbf{A}_i = \mathbf{LDL}' \)$

So,

\[ \mathbf{A}_i^{-1} = (\mathbf{L}')^{-1}\mathbf{D}^{-1}\mathbf{L}^{-1} \]
\[\begin{split} \begin{align} (\mathbf{L}')^{-1}\mathbf{D}^{-1}\mathbf{L}^{-1} &= \begin{bmatrix} \mathbf{I} & -\mathbf{q}_i \\ \mathbf{0} & 1 \end{bmatrix} \begin{bmatrix} \mathbf{A}_{i-1}^{-1} & \mathbf{0} \\ \mathbf{0} & (a_{ii} - \mathbf{q}'_i\mathbf{A}_{i-1}\mathbf{q}_i)^{-1} \end{bmatrix} \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ -\mathbf{q}'_i & 1 \end{bmatrix} \\ &= \begin{bmatrix} \mathbf{A}_{i-1}^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} + \begin{bmatrix} -\mathbf{q}_i\\ 1 \end{bmatrix} a^{ii} \begin{bmatrix} -\mathbf{q}'_i & 1 \end{bmatrix} \end{align} \end{split}\]

where \(a^{ii}= (a_{ii} - \mathbf{q}'_i\mathbf{A}_{i-1}\mathbf{q}_i)^{-1}\). As shown below, this can be computed without computing \(\mathbf{A}_{i-1}\), which can be a large dense matrix.

4.7. Calculation of \(a^{ii}\)#

When both parents are known, \(\mathbf{q}'_i\mathbf{A}_{i-1}\mathbf{q}_i = \frac{1}{4}(a_{m_im_i} + a_{f_if_i} + 2a_{m_if_i})\), and because \(a_{ii} = 1 + \frac{1}{2}a_{m_if_i}\),

\[\begin{split} \begin{align} a^{ii} &= [1 - \frac{1}{4}(a_{m_im_i} + a_{f_if_i})]^{-1} \\ &= (\frac{1}{2} - \frac{1}{4}f_m - \frac{1}{4}f_f)^{-1}\\ &= [\frac{1}{4}(2 - f_m - f_f)]^{-1}\\ &= \frac{4}{2 - f_m - f_f}, \end{align} \end{split}\]

where \(f_m\) and \(f_f\) are the inbreeding values of the parents \(m\) and \(f\) of \(i\).

4.8. Algorithm#

  • Set \(\mathbf{A}^{-1} = \mathbf{0}\)

  • Compute \(a^{ii}\) for all animals

  • For each animal add to \(\mathbf{A}^{-1}\)

    • \(a^{ii}\) to (i,i)

    • -1/2\(a^{ii}\) to \((i,f_i)\), \((f_i,i)\),\((i,m_i)\),\((m_i,i)\)

    • 1/4\(a^{ii}\) to \((f_i,f_i)\), \((f_i,m_i)\),\((m_i,f_i)\),\((m_i,m_i)\)

4.9. Numerical Example#

using LinearAlgebra
A2i = I(2)*1.0

2×2 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅ 
  ⋅   1.0

q = [0.5, 0.5]
a33i = 1/(A[3,3] - q'A2*q)
A3i  = [A2i zero(q)
        zero(q')   0] + [-q; 1]*[-q' 1].*a33i

3×3 Matrix{Float64}:
  1.5   0.5  -1.0
  0.5   1.5  -1.0
 -1.0  -1.0   2.0

inv(A[1:3,1:3])

3×3 Matrix{Float64}:
  1.5   0.5  -1.0
  0.5   1.5  -1.0
 -1.0  -1.0   2.0

q    = [0.5, 0.5, 0]
a44i = 1/(A[4,4] - q'A3*q)
A4i  = [A3i zero(q)
        zero(q')   0] + [-q; 1]*[-q' 1].*a44i

4×4 Matrix{Float64}:
  2.0   1.0  -1.0  -1.0
  1.0   2.0  -1.0  -1.0
 -1.0  -1.0   2.0   0.0
 -1.0  -1.0   0.0   2.0

inv(A[1:4,1:4])

4×4 Matrix{Float64}:
  2.0   1.0  -1.0  -1.0
  1.0   2.0  -1.0  -1.0
 -1.0  -1.0   2.0  -0.0
 -1.0  -1.0   0.0   2.0

q    = [0, 0, 0.5, 0.5]
[-q; 1]*[-q' 1]

5×5 Matrix{Float64}:
  0.0   0.0   0.0    0.0   -0.0
  0.0   0.0   0.0    0.0   -0.0
  0.0   0.0   0.25   0.25  -0.5
  0.0   0.0   0.25   0.25  -0.5
 -0.0  -0.0  -0.5   -0.5    1.0

a55i = 1/(A[5,5] - q'A4*q)
A5i  = [A4i zero(q)
        zero(q')   0] + [-q; 1]*[-q' 1].*a55i

5×5 Matrix{Float64}:
  2.0   1.0  -1.0  -1.0   0.0
  1.0   2.0  -1.0  -1.0   0.0
 -1.0  -1.0   2.5   0.5  -1.0
 -1.0  -1.0   0.5   2.5  -1.0
  0.0   0.0  -1.0  -1.0   2.0

inv(A)

5×5 Matrix{Float64}:
  2.0   1.0  -1.0  -1.0  -0.0
  1.0   2.0  -1.0  -1.0  -0.0
 -1.0  -1.0   2.5   0.5  -1.0
 -1.0  -1.0   0.5   2.5  -1.0
  0.0   0.0  -1.0  -1.0   2.0

4.10. Inverse of any positive definite symmetric \(\mathbf{A}\)#

\(\mathbf{A}\) is an \(n\times n\) positive definite matrix that is partitioned as

\[\begin{split} \begin{align*} \begin{bmatrix} \mathbf{A}_{11} & \mathbf{a}_{12} \\ \mathbf{a}_{21} & \mathbf{a}_{22} \end{bmatrix} &= \begin{bmatrix} \mathbf{A}_{11} & \mathbf{A}_{11}\mathbf{q} \\ \mathbf{q}'\mathbf{A}_{11} & \mathbf{q}'\mathbf{A}_{11}\mathbf{q} + \mathbf{r} \end{bmatrix}\\ \end{align*} \end{split}\]

where \(\mathbf{A}_{11}\) is \(n-1\times n-1\), \(\mathbf{a}_{12}\) is \(n-1\times 1\) and \(\mathbf{a}_{22}\) is \(1\times 1\).

So,

\[ \mathbf{q} = \mathbf{A}_{11}^{-1}\mathbf{a}_{12}, \]
\[ r = a_{22} - \mathbf{q}'\mathbf{A}_{11}\mathbf{q} \]

and

  • if \(n=1\)

\[ \mathbf{A}^{-1} = \frac{1}{\mathbf{A}} \]

  • else:

\[\begin{split} \mathbf{A}^{-1} = \begin{bmatrix} \mathbf{A}_{11}^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} + \begin{bmatrix} \mathbf{q}\mathbf{q}' & -\mathbf{q} \\ -\mathbf{q}' & 1 \end{bmatrix}\frac{1}{r} \end{split}\]

4.11. Recursive function for inverse#

function invA(A)
    n  = size(A,1)
    if n==1
        return 1/A[1,1]
    else
        Ai = zeros(n,n) # inverse of A
        s1 = 1:n-1
        s2 = n
        A11 = A[s1,s1]
        A11i = invA(A11)
        A12 = A[s1,s2]
        A22 = A[s2,s2]
        q = A11i*A12
        r = A22 - q'A11*q
        ri = 1/r
        Ai = [A11i .+ ri*q*q'    -q*ri
              -q'*ri              ri]
        return Ai
   end
end

invA (generic function with 1 method)

4.12. Numerical Example#

A = 5
invA(A)

0.2

X = randn(10,5)
A = X'X
Ai = invA(A)
round.(A*Ai,digits=3)

5×5 Matrix{Float64}:
  1.0  -0.0  -0.0  -0.0  -0.0
 -0.0   1.0  -0.0  -0.0  -0.0
  0.0   0.0   1.0   0.0   0.0
 -0.0   0.0  -0.0   1.0   0.0
 -0.0  -0.0   0.0  -0.0   1.0