随着基因分型技术的发展，成本随之降低^[1]，基因组选择逐渐成为畜禽育种的主流方法^[2-9]。当前基因组选择应用最多的评估方法是一步基因组最佳线性无偏预测(ssGBLUP)^[10-16]，该方法在基因组最佳线性无偏预测(genome best linear unbiased prediction, GBLUP)方法的基础上将基因组分型和未分型的个体统一进行遗传评估^[17-19]，能进一步提高遗传评估的准确性^[20-21]，且和常规育种软件无缝衔接^[22-24]。GBLUP和ssGBLUP都涉及到基因组亲缘矩阵(G)的求逆运算^[25-26]。大规模矩阵求逆是一个十分耗时的过程^[27]，而育种具有时效性，随着数据的积累，快速进行基因组估计育种值(GEBV)估计是基因组选择亟待解决的问题。

大型稠密矩阵求逆对于计算设备的要求很高，计算时间长。在寄希望于计算设备发展更快、成本更低的同时，数量遗传学家也在尝试寻找更加高效、快速的算法，来实现G矩阵的求逆。更新法^[28](updating the inverse)、APY法^[29-30](algorithm for proven and young animals)、迭代法^[31](recursive algorithm for approximation of inverse)等一系列方法都在一定程度上实现了G矩阵的快速求逆。更新法将新的基因分型个体增加到G^-1中，避免了重新构建G矩阵与求逆的过程。该方法在实际生产中具有一定的应用价值，但并没有解决求逆效率问题。APY法将基因分型后的个体分为核心个体和非核心个体，利用分块矩阵的思想^[32]将非核心个体的子块转化为一个对角矩阵，在损失一定准确性基础上提高了求逆的效率。张国亮等^[33]提出了大型实对称矩阵分块迭代求逆算法，每次迭代一行元素，避免了矩阵直接求逆的过程，但大量的矩阵、向量相乘影响了大型矩阵的求逆效率。

本研究旨在在分块迭代求逆算法基础上寻找一种更高效的求逆方法，从而有效提高GEBV的估计效率。

1 材料与方法 1.1 数据来源与计算设备

利用VanRaden^[34]报道的亲缘矩阵构建方法计算G矩阵，并通过与分子亲缘矩阵加权求和来避免矩阵奇异，G^*=wG+(1-w)A，G^*为基因组亲缘矩阵和分子亲缘矩阵加权求和后的矩阵，w为权重，其值可取0.95或0.98^[35]；A为分子亲缘矩阵。本研究利用该方法对基因组亲缘矩阵及其对角子块进行加权处理，保证矩阵和各子块可逆，以便于后续方法验证。基因型数据分为两部分，其中4 000阶、8 000阶、12 000阶、16 000阶、32 000阶、64 000阶和128 000阶的矩阵所用基因型数据是使用QMSIM^[36]模拟软件产生的模拟数据；剩余数据为1 734头猪的Geneseek 50K芯片数据，来源于各个猪场(四川江油新希望(971头)、浙江天蓬集团(604)、四川绵阳明兴农业(159头))。计算设备分为两部分：一是自用台式电脑，CPU为intel i7-8550U 1.8 GHz，内存8 GB；二是学校的公共服务器，CPU为intel Xeon E5-2630 2.6 GHz。

1.2 方法

1.2.1 方法1 　　分块迭代方法可以避免直接对矩阵进行求逆运算，但是其迭代的次数太多，每次迭代都涉及到矩阵的分块和组合以及子块间的乘法。减少迭代次数理论上能够提高求逆效率，同时引入intel数学核心函数库(MKL)解决子块直接求逆问题并进一步提高求逆效率。每一次分块都将p_t分为一个2×2、3×3或者更大的矩阵(方法1)，这样能有效减少迭代次数。公式：

$ \begin{array}{c} \mathit{\boldsymbol{W}}_t^{ - 1} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{W}}_{t - 1}}}&{{\mathit{\boldsymbol{r}}_t}}\\ {\mathit{\boldsymbol{r}}_t^T}&{{\mathit{\boldsymbol{p}}_t}} \end{array}} \right]^{ - 1}} = \\ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}_{t - 1}^{ - 1} + \mathit{\boldsymbol{W}}_{t - 1}^{ - 1}{\mathit{\boldsymbol{r}}_t}\mathit{\boldsymbol{\beta }}_t^{ - 1}\mathit{\boldsymbol{r}}_t^T\mathit{\boldsymbol{W}}_{t - 1}^{ - 1}}&{ - \mathit{\boldsymbol{W}}_{t - 1}^{ - 1}{\mathit{\boldsymbol{r}}_t}\mathit{\boldsymbol{\beta }}_t^{ - 1}}\\ { - \mathit{\boldsymbol{\beta }}_t^{ - 1}\mathit{\boldsymbol{r}}_t^T{\mathit{\boldsymbol{W}}_{t - 1}}}&{\mathit{\boldsymbol{\beta }}_t^{ - 1}} \end{array}} \right] \end{array} $

(1)

$ {\mathit{\boldsymbol{\beta }}_t} = {\mathit{\boldsymbol{p}}_t} - \mathit{\boldsymbol{r}}_t^T\mathit{\boldsymbol{W}}_{t - 1}^{ - 1}{\mathit{\boldsymbol{r}}_t} $

(2)

W_t^-1为第t次迭代求逆完成后的逆矩阵，t表示迭代的次数，T表示转置矩阵。W_t-1、r_t、p_t和r_t^T为第t次迭代所用子块，β_t为第t次迭代后逆矩阵的右下角子块。

1.2.2 方法2 　　对于一个N (N＞0)阶矩阵，方法1中的p_t最大阶数为N/2，此时矩阵被分为4个大小相等的矩阵，重复对左上角子块进行分块降阶，理论上可以进一步提高矩阵求逆的效率。进行m (m＞0)次分块后需要直接求逆的矩阵的阶数为$ \frac{N}{{{2^m}}}$。

$ {\mathit{\boldsymbol{G}}^{* - 1}} = {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{a}}&\mathit{\boldsymbol{b}}\\ \mathit{\boldsymbol{c}}&\mathit{\boldsymbol{d}} \end{array}} \right]^{ - 1}} $

(3)

$ {\mathit{\boldsymbol{a}}^{ - 1}} = {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{p}}&\mathit{\boldsymbol{q}}\\ \mathit{\boldsymbol{k}}&\mathit{\boldsymbol{l}} \end{array}} \right]^{ - 1}} $

(4)

G^*-1为加权基因组亲缘矩阵的逆矩阵，a、b、c、d为G^*矩阵的子块，a^-1为子块a的逆矩阵，p、q、k、l为子块a的子矩阵块。两次分块后矩阵p的阶数为N/4。然后根据公式(1)和(2)求得a^-1，继而求得G^*-1。

2 结果 2.1 方法1结果

根据张国亮等^[33]的方法，用真实数据测试矩阵求逆的效率，可以看到随着子块p_t阶数增加，矩阵求逆的效率明显提升。详细结果见图 1 (计算设备：台式电脑)。

2.2 方法2结果

将矩阵分成4个大小相等的子块，再对右上角的子块继续分块或直接求逆。为了使结果对比更加明显，本研究使用模拟数据来测试，并分开画图。可以发现，随着矩阵阶数增加，方法2的优势愈加明显，详细结果见图 2、图 3(计算设备：服务器)。

方法2在计算矩阵求逆上有更高的效率，4 000× 4 000的基因组亲缘矩阵求逆效率是inv方法的1.084倍；128 000×128 000的矩阵，与inv方法的加速比为1.805倍。但是计算效率并不会随迭代次数无限增加，64 000阶矩阵求逆效率随迭代次数先升高后降低，详细结果见图 4(计算设备：服务器)。

2.3 方法1与方法2结果对比

方法1在计算4 000×4 000的基因组亲缘矩阵逆矩阵时，与inv方法的加速比为0.898。而16 000× 16 000的计算速度为inv方法的1.006倍。当方法1中p_t的阶数等于方法2子块的阶数时，两种方法的计算思路和计算步骤完全一致，理论上两种方法在计算时间上应该非常接近。在实际测试中两种方法的计算效率存在差异，详细结果如图 5 (计算设备：台式电脑)。

3 讨论

利用公式G^*=wG+(1-w)A对G矩阵处理后，可以避免矩阵奇异。当子块不可逆时，同样利用该公式对子块进行处理，A为分子亲缘矩阵中对应子块个体的元素，该方式能有效避免子块奇异且不影响后续评估GEBV^[35]。

在张国亮等^[33]的报道中，700×700的矩阵求逆时间为3.506 s (计算设备：intel i5-3470)，而在本研究中，方法1(n=3)计算885×885大小矩阵的逆时间为0.365 s。Meyer等^[28]对10 000×10 000的矩阵直接求逆，消耗45 s，再用更新法在前面的逆矩阵基础上增加2 000阶，消耗44 s。最终完成一个12 000×12 000大小的矩阵求逆，消耗时间89 s(计算设备：intel i7-3930K)。在本研究中，方法2(m=3)计算12 000×12 000的矩阵逆消耗时间10.149 s。本研究的方法计算效率更高，但计算设备不同，因此该部分仅供参考。

通过图 2和图 3来看，方法2比直接调用LAPACKE_dgetri函数具有一定的优势，且随着迭代次数的增加优势更加明显。对于一个N阶矩阵，矩阵乘法和矩阵求逆的复杂度都为O(N³)。基因组亲缘矩阵为对称矩阵，结合公式(1)和公式(2)，一次分块后，方法2需要计算两次子块的求逆、4次子块间的乘法。因此，方法2(m=1)的复杂度为$ O\left( {{{\left( {\frac{N}{2}} \right)}^3} \times 6} \right) = O\left( {\frac{{3{N^3}}}{4}} \right)$。所以，理论上方法2比直接求逆效率更高，实际结果符合预期。

方法2要求被分块的矩阵为偶数，这在实际应用中存在一定限制。方法1可以任意分块，更符合实际情况。通过图 5可以看出，选择合适的p_t阶数，方法1也能超越MKL函数直接求逆的效率，但和方法2 (m=1)相比仍然较低。两种方法在思路和步骤上几乎完全一致，但在分块方式上存在一定差异。方法1首先确定了子块d (公式(3))的大小，子块b的大小会随迭代逐渐增大。为了保证可以任意分块，子块d和b的赋值过程必须单独进行。方法2每次分块后子块d和b的阶数相同，因此赋值过程可以在同一个循环中完成。相比方法2，方法1每次迭代过程都要多出一个循环赋值的过程，这是方法2比方法1效率更高的原因。但方法2随着分块次数增加，矩阵求逆的效率先增加后降低(图 4)。可能原因是随着分块次数的增加，矩阵的分块与子块的组合消耗了过多的时间。对于不同大小的矩阵，其最佳分块次数是不同的，64 000×64 000的矩阵最佳分块次数是四次(图 4)，但128 000×128 000 (图 3)的矩阵在第二次分块后效率开始降低。从图 2、图 3和图 4来看，最佳的分块次数出现在2~4次之间，最佳分块次数及其原理是下一步研究的方向。

4 结论

综上所述，方法1与方法2在效率上都能超越MKL中的直接求逆函数。方法2在效率上更具有优势，但其分块方式在实际应用中存在限制。在实际应用中，两种方法可以选择性使用。对于两种方法的进一步优化，应该从数据结构方面进行，降低内存的消耗、减少数据赋值的次数，进一步提高矩阵求逆的效率。