在 ssGBLUP 中需要构建 A22 矩阵，通过 Colleau 算法可以避免直接构建 A 矩阵，节约内存和时间。

A 阵分解

我们可以将 \(\mathbf{A}\) 分解为 \(\mathbf{A} = \mathbf{(I-P)^{-1}D(I-P)^{-1\prime}}\) （要求系谱按照出生日期顺序排序）。

其中 \(\mathbf{D}\) 是一个对角矩阵，其元素计算公式如下，其中 \(F_{i}\) 是近交系数。 \[ \begin{aligned} D_{ii} & = 0.5 - 0.25(F_{s_{i}} + F_{d_{i}}) \quad \text{当 i 的父母均已知}\\ & = 0.75 - 0.25 F_{p_{i}} \quad \text{当 i 的单个亲本已知}\\ & = 1 \quad \text{当 i 的两个亲本均未知}\\ \end{aligned} \] 这里 \(\mathbf{P}\) 是一个下三角矩阵，对角线元素为 0，非对角线只有亲子间为 0.5，其余元素也为 0 。

\(\mathbf{I-P}\) 矩阵举例如下： \[ \left[\begin{array}{cccccc} 1 & & & & & \\ 0 & 1 & & & \\ -0.5 & 0 & 1 & & \\ -0.5 & -0.5 & 0 & 1 & \\ 0 & 0 & -0.5 & -0.5 & 1 & \\ -0.5 & 0 & 0 & -0.5 & 0 & 1 \\ 0 & 0 & 0 & 0 & -0.5 & -0.5 & 1 \end{array}\right] \]

Colleau 算法

Colleau (2002) 年证明 \(\mathbf{A}\) 阵与一个向量 \(\mathbf{t}\) 的乘积可以在一个线性的时间内计算完成，如下 \[ \mathbf{v} = \mathbf{At} = \mathbf{(I-P)^{-1}D[(I-P)^{-1\prime}t]} \] 这里我们可以分步骤进行计算：

• 首先我们定义 \(\mathbf{r} = \mathbf{(I-P)^{-1\prime}t}\) ，转换一下我们可以通过求解三角方程组 \(\mathbf{(I-P)^{\prime}r=t}\) 得到 \(\mathbf{r}\) 。

• 然后我们得到 \(\mathbf{v} = \mathbf{(I-P)^{-1}Dr}\) ，转换一下我们可以通过求解三角方程组 \(\mathbf{(I-P)v=Dr}\) 得到 \(\mathbf{v}\) 。

标量公式如下（第一行是采用 outer-product 的原理；第二行是采用向量内积的原理） \[ r_{i} = r_{i}+t_{i}; \quad r_{si} = r_{si} + r_{i}/2; \quad r_{di} = r_{di} + r_{i}/2; \quad i=n,\cdots,1 \]

\[ v_{i} = d_{i}r_{i}+ (v_{si}+v_{di})/2, \quad i = 1,\cdots,n \]

这里的 \(si\) 和 \(di\) 表示个体 \(i\) 的父母本。

通过 Colleau 算法可以用于计算子矩阵，举个例子，下面展示了如何计算 \(\mathbf{A}_{11}\mathbf{q}\)，\(\mathbf{A}_{21}\mathbf{q}\) 。 \[ \left[\begin{array}{c} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{21} & \mathbf{A}_{22} \\ \end{array}\right] \left[\begin{array}{c} \mathbf{q} \\ \mathbf{0} \\ \end{array}\right] = \left[\begin{array}{c} \mathbf{A}_{11} \mathbf{q} \\ \mathbf{A}_{21} \mathbf{q} \\ \end{array}\right] \]

构建 A22 矩阵

这里需要注意的是，构建 \(\mathbf{A}\) 阵时是按照出生日期顺序排序，但是我们提取A22 矩阵要求按照无基因型个体在前，有基因型个体在后的顺序排列，因此我们还要对原始的 \(\mathbf{A}\) 阵的行和列进行重新重排。我们设相应的置换矩阵为 \(\mathbf{Q}\) ，重排后按照无基因型个体在前，有基因型个体在后的顺序的矩阵称为 \(\mathbf{A}^{*}\) 矩阵，则有 \[ \mathbf{A}^{*} = \mathbf{QAQ^{\prime}} = \mathbf{Q(I-P)^{-1}D(I-P)^{-1\prime}Q^{\prime}} \] 此时我们逐列提取 A22 矩阵的每一列，举个例子，假如基因型个体在 \(\mathbf{A}^{*}\) 矩阵中的位置是 80 到 100，那么我们构建第一个 \(\mathbf{t}\) 向量就是在 80位置为1，其它位置为0；构建第二个 \(\mathbf{t}\) 向量就是在 81位置为1，其它位置为0 …以此类推。

通过 \(\mathbf{A}^{*}\mathbf{t}\) 便可以得到 A22 矩阵个体的一列，其计算公式为 \[ \mathbf{v}^{*} = \mathbf{A}^{*}\mathbf{t} = \mathbf{Q(I-P)^{-1}D(I-P)^{-1\prime}Q^{\prime}t} \] 这里我们可以同样分步骤进行计算：

• 首先计算 \(\mathbf{t}^{*} = \mathbf{Q^{\prime}t}\)

• 我们定义 \(\mathbf{r} = \mathbf{(I-P)^{-1\prime}t^{*}}\) ，转换一下我们可以通过求解三角方程组 \(\mathbf{(I-P)^{\prime}r=t^{*}}\) 得到 \(\mathbf{r}\) 。

• 然后我们得到 \(\mathbf{v} = \mathbf{(I-P)^{-1}Dr}\) ，转换一下我们可以通过求解三角方程组 \(\mathbf{(I-P)v=Dr}\) 得到 \(\mathbf{v}\) 。

• 然后我们通过 \(\mathbf{v}^{*} = \mathbf{Qv}\) 得到 \(\mathbf{v}^{*}\) 。

最后这里得到的 \(\mathbf{v}^{*}\) 还包含一些无关元素，按照上例，还需要从 \(\mathbf{v}^{*}\) 提取从 80 到 100 的位置的元素，这才是真正的 A22 矩阵的一列元素。

参考文献

1. Colleau J J. An indirect approach to the extensive calculation of relationship coefficients[J]. Genetics Selection Evolution, 2002, 34(4): 409.

2. Misztal I, Aguilar I, Legarra A, et al. A unified approach to utilize phenotypic, full pedigree, andgenomic information for genetic evaluation[C]//9. World Congress on Genetics Applied to Livestock Production. 2010.