当数据量很大时，对混合模型方程组的系数矩阵直接求逆困难或不可实现，需要采用迭代求解的方法。

迭代求解的基本方法

雅可比迭代法

设有线性方程组 \[ \mathbf{Ab} = \mathbf{r} \] 其中，\(\mathbf{A}\) 是 \(p \times p\) 的系数矩阵；\(\mathbf{b}\) 为待求解的 \(p\) 维向量；\(\mathbf{r}\) 为右手项的 \(p\) 维向量。

如果见矩阵 \(\mathbf{A}\) 分解为 \[ \mathbf{A} = \mathbf{M+N} \] 其中，\(\mathbf{M}\) 为一可逆矩阵，将其代入上式，则有 \[ \mathbf{Mb} = \mathbf{-Nb+r} \] 或写成 \[ \mathbf{b} = \mathbf{Qb +f} \] 其中，\(\mathbf{Q} = \mathbf{-M^{-1}N}\)，\(\mathbf{f} = \mathbf{M^{-1}r}\) 。

那么我们可以构建出迭代方法 \[ \mathbf{b^{(t+1)}} = \mathbf{Qb^{(t)}+f} \] 如果在 \(t \rightarrow \infty\) 时 , \(\mathbf{b^{(t)}} \rightarrow \mathbf{b^{*}}\) ，我们称迭代收敛，且 \(\mathbf{b^{*}}\) 就是该方程组的解。

例如，对于 \(p =3\) 时，设 \(\mathbf{A} = \{a_{ij}\}\) ，其对角线元素均不为0，即 \(a_{ii} \neq 0\) ，则我们可取 \(\mathbf{M}\) 和 \(\mathbf{N}\) 为 \[ \mathbf{M}=\left[\begin{array}{ccc} a_{11} & 0 & 0 \\ 0 & a_{22} & 0 \\ 0 & 0 & a_{33} \end{array}\right], \mathbf{N}=\left[\begin{array}{ccc} 0 & a_{12} & a_{13} \\ a_{21} & 0 & a_{23} \\ a_{31} & a_{32} & 0 \end{array}\right] \] 那么我们可以得到 $$ \[\begin{aligned} \mathbf{Q} &= \mathbf{-M^{-1}N} \\ &=- \left[\begin{array}{ccc} 1/a_{11} & 0 & 0 \\ 0 & 1/a_{22} & 0 \\ 0 & 0 & 1/a_{33} \end{array}\right] \left[\begin{array}{ccc} 0 & a_{12} & a_{13} \\ a_{21} & 0 & a_{23} \\ a_{31} & a_{32} & 0 \end{array}\right] \\ &= \left[\begin{array}{ccc} 0 & -\frac{a_{12}}{a_{11}} & -\frac{a_{13}}{a_{11}} \\ -\frac{a_{21}}{a_{22}} & 0 & -\frac{a_{23}}{a_{22}} \\ -\frac{a_{31}}{a_{33}} & -\frac{a_{32}}{a_{33}} & 0 \end{array}\right] \\ \\ \mathbf{f} &= \mathbf{M^{-1}r} \\ &= \left[\begin{array}{ccc} 1/a_{11} & 0 & 0 \\ 0 & 1/a_{22} & 0 \\ 0 & 0 & 1/a_{33} \end{array}\right] \left[\begin{array}{ccc} r_{1} \\ r_{2} \\ r_{3} \end{array}\right] \\ &= \left[\begin{array}{l} \frac{r_{1}}{a_{11}} \\ \frac{r_{2}}{a_{22}} \\ \frac{r_{3}}{a_{33}} \end{array}\right] \\ \end{aligned}\] \[ 因此，我们有方程组 \] =+\[ 即 \] { \[\begin{array}{l} b_{1}=\left(r_{1}-a_{12} b_{2}-a_{13} b_{3}\right) / a_{11} \\ b_{2}=\left(r_{2}-a_{21} b_{1}-a_{23} b_{3}\right) / a_{22} \\ b_{3}=\left(r_{3}-a_{31} b_{1}-a_{32} b_{2}\right) / a_{33} \end{array}\] . \[ 因此，我们得到了一个迭代公式 \] { \[\begin{array}{l} b_{1}^{(t+1)}=\left(r_{1}-a_{12} b_{2}^{(t)}-a_{13} b_{3}^{(t)}\right) / a_{11} \\ b_{2}^{(t+1)}=\left(r_{2}-a_{21} b_{1}^{(t)}-a_{23} b_{3}^{(t)}\right) / a_{22} \\ b_{3}^{(t+1)}=\left(r_{3}-a_{31} b_{1}^{(t)}-a_{32} b_{2}^{(t)}\right) / a_{33} \end{array}\]

. \[ 将上面的情形推广至一般情况，即可以将 $\mathbf{Ab} = \mathbf{r}$ 改写为 \] b_{i}=(r_{i}-{j=1}^{i-1} a{i j} b_{j}-{j=i+1}^{p} a{i j} b_{j}) / a_{i i}, (i=1,2,,p) \[ 由此可得如下的一般迭代式 \] b_{i}^{{(t+1)}=(r_{i}-_{j=1}}{i-1} a_{i j} b_{j}^{{(t)}-_{j=i+1}}{p} a_{i j} b_{j}^{(t)}) / a_{i i}, (i=1,2,,p) \[ 或者写成 \] b_{i}^{{(t+1)}=b_{i}}{(t)}+(r_{i}-{j=1}^{p} a{i j} b_{j}^{(t)}) / a_{i i}, i=1, , p \[ 写成矩阵的形式，我们有 \] = \[ 或者 \] = $$ 其中 \(\mathbf{D} = Diag\{\mathbf{A}\}\) ，是由 \(\mathbf{A}\) 中的对角线元素构成的对角线矩阵。

这个迭代法称为 雅可比 (Jacobi) 迭代法。在计算时一般按 \(b_{1}^{(t+1)}, b_{2}^{(t+1)}, \cdots\), \(b_{p}^{(t+1)}\) 的顺序进行。注意在计算 \(b_{i}^{(t+1)}\) 时, 在它前面的 \(b_{1}^{(t+1)}, \cdots, b_{i-1}^{(t+1)}\) 已经计算出 来了, 而雅可比法并不利用这些最新的迭代值, 仍用 \(b_{1}^{(t)}, \cdots, b_{i-1}^{(t)}\) 。

写成伪代码的格式为 \[ \begin{aligned} \text{do i=1,n}& \\ &b_{i}^{(t+1)}=b_{i}^{(t)}+\left(r_{i}-\sum_{j=1}^{p} a_{i j} b_{j}^{(t)}\right) / a_{i i} \\ \text{end do}&\\ \end{aligned} \]

另一种理解

假设系数矩阵的对角线元素均不为0 （如果存在某些行的对角线元素为 0 ，可以与邻近的行进行交换位置），那么雅可比迭代法其实就是用第一个方程计算 \(\hat{b}_{1}\) ，用第二个方程计算 \(\hat{b}_{2}\) ，依次类推。以上面的 \(p=3\) 的例子为例，我们将方程组写成标量形式，即： \[ \begin{aligned} &a_{11} b_{1}+a_{12} b_{2}+a_{13} b_{3}=r_{1} \\ &a_{21} b_{1}+a_{22} b_{2}+a_{23} b_{3}=r_{2} \\ &a_{31} b_{1}+a_{32} b_{2}+a_{33} b_{3}=r_{3} \end{aligned} \] 用第一个方程计算 \(\hat{b}_{1}\) ，用第二个方程计算 \(\hat{b}_{2}\) ，用第三个方程计算 \(\hat{b}_{3}\) ，得到下式 \[ \begin{aligned} &b_{1}^{(t+1)}=\left(1 / a_{11}\right)\left(r_{1}-a_{12} b_{2}^{(t)}-a_{13} b_{3}^{(t)}\right) \\ &b_{2}^{(t+1)}=\left(1 / a_{22}\right)\left(r_{2}-a_{21} b_{1}^{(t)}-a_{23} b_{3}^{(t)}\right) \\ &b_{3}^{(t+1)}=\left(1 / c_{33}\right)\left(r_{3}-c_{31} b_{1}^{(t)}-c_{32} b_{2}^{(t)}\right) \end{aligned} \] 我们同样得到雅可比的迭代公式，这种方式更好理解一点。

GS迭代法

如果我们对这一点进行修改, 即在每个 \(b_{i}^{(u+1)}\) 的计算中尽量利用最新的迭代值, 于是有 \[ b_{i}^{(t+1)}=\left(r_{i}-\sum_{j=1}^{i-1} a_{i j} b_{j}^{(t+1)}-\sum_{j=i+1}^{p} a_{i j} b_{j}^{(t)}\right) / a_{i i} \] 或者 \[ b_{i}^{(t+1)}=b_{i}^{(t)}+\left(r_{i}-\sum_{j=1}^{i-1} a_{i j} b_{j}^{(t+1)}-\sum_{j=i}^{p} a_{i j} b_{j}^{(t)}\right) / a_{i i} \] 这种迭代方法称为 高斯-塞德尔 (Gauss-Seidel, GS) 迭代法。此时编程更加容易，因为新解可以直接覆盖旧的解。伪代码如下 \[ \begin{aligned} \text{do i=1,n}& \\ &b_{i}=b_{i}+\left(r_{i}-\sum_{j=1}^{p} a_{i j} b_{j}\right) / a_{i i} \\ \text{end do}&\\ \end{aligned} \] 上式可以改写为 $$ \[\begin{aligned} b_{i}^{(t+1)}&=\left(r_{i}-\sum_{j=1}^{i-1} a_{i j} b_{j}^{(t+1)}-\sum_{j=i+1}^{p} a_{i j} b_{j}^{(t)}\right) / a_{i i} \\ a_{i i} b_{i}^{(t+1)} &=\left(r_{i}-\sum_{j=1}^{i-1} a_{i j} b_{j}^{(t+1)}-\sum_{j=i+1}^{p} a_{i j} b_{j}^{(t)}\right) \\ \sum_{j=1}^{i} a_{i j} b_{j}^{(t+1)} &= r_{i}-\sum_{j=i+1}^{p} a_{i j} b_{j}^{(t)} \\ \end{aligned}\]

\[ 我们看如何将其写成矩阵的形式，首先我们矩阵 $\mathbf{A}$ 分解为如下形式 \] = $$ 其中 \(\mathbf{D} = Diag\{\mathbf{A}\}\) ，是由 \(\mathbf{A}\) 中的对角线元素构成的对角线矩阵。 \(\mathbf{L, U}\) 均为严格三角矩阵: \(\mathbf{L}\) 为严格下 三角矩阵，即对角线以下的元素； \(\mathbf{U}\) 为严格上三角矩阵，即对角线以上的元素。

于是我们有 \[ \begin{aligned} \mathbf{(L+D) \mathbf{b^{(t+1)}}} &= \mathbf{r-U b^{(t)}} \\ \mathbf{b^{(t+1)}} &= \mathbf{(L+D)^{-1}\left\{r-U b^{(t)}\right\}} \\ &= \mathbf{(L+D)^{-1}\left\{r-[A-(L+D)] b^{(t)}\right\}}\\ \end{aligned} \] 与雅可比迭代的式子比较发现，这里是用 \(\mathbf{L+D}\) 替代了雅可比迭代式子中的 \(\mathbf{D}\) ，由于 \(\mathbf{L+D}\) 比 \(\mathbf{D}\) 更接近 \(\mathbf{A}\) , 因此 GS 迭代法的收敛速度一般更快。

但要注意的是，这两种方法并不是在任何情况下都能收敛的。对于 GS 迭代法，它能够收敛的一个充分条件是方程组的系数矩阵为正定矩阵或半正定矩阵，由于混合线性方程组的系数矩阵一般来说都是正定矩阵，所以基本上可以不用担心它们的收敛问题。对于雅可比迭代法，一般来说，对于系数矩阵为对角优势矩阵（即在每一行中，均有对角线元素的值大于所有非对角线元素的值）的方程组，它是可以收敛的。但是混合模型方程组中的系数矩阵往往不是对角优势矩阵，因此不能保证它收敛，因此只有在对数据文件迭代 (iteration on data) 中才可能用到雅可比迭代。

SOR 迭代法

GS 迭代法有时收敛速度较慢，此时可考虑逐次超松弛迭代法 (successive over-relaxation, SOR) ，其迭代公式为 \[ b_{i}^{(t+1)}=b_{i}^{(t)}+ \omega \left(r_{i}-\sum_{j=1}^{i-1} a_{i j} b_{j}^{(t+1)}-\sum_{j=i}^{p} a_{i j} b_{j}^{(t)}\right) / a_{i i} \] 其中，\(\omega\) 是松弛因子，适当选取 \(\omega\) 可加快收敛速度，通常取 \(1<\omega<2\) （那么如何选择一个合适的 \(\omega\) 就成了一个新的问题）。当 \(\omega = 1\) 时，SOR 迭代法就变成 GS 迭代法。 SOR 方法的伪代码如下 \[ \begin{aligned} \text{do i=1,n}& \\ &b_{i}=b_{i}+\omega\left(r_{i}-\sum_{j=1}^{p} a_{i j} b_{j}\right) / a_{i i} \\ \text{end do}&\\ \end{aligned} \]

PCG

在复杂问题中，预处理共轭梯度法 (Preconditioned conjugate gradient, PCG) 比前面三种方法收敛速度都要快，其理论可见 Meurant (1999) 。

PCG 和二阶雅可比算法有一些相似，二阶雅可比算法如下 (在动物模型，合适的参数是 \(\alpha = 0.7-0.9\)， \(\beta =1\) 。) \[ \mathbf{b}^{(\mathrm{n}+1)}=\mathbf{b}^{(\mathrm{n})}+\alpha\left(\mathbf{b}^{(\mathrm{n})}-\mathbf{b}^{(\mathrm{n}-1)}\right)+\beta \mathbf{D}^{-1}\left(\mathbf{r}-\mathrm{A} \mathbf{b}^{(\mathrm{n})}\right) \] PCG 算法采用相似的公式，只是 \(\alpha\) 和 \(\beta\) 在每一轮迭代中需要重新计算，并且这里的 \(\mathbf{D}\) 改为了 \(\mathbf{M}\) 矩阵，称为预处理矩阵。\(\mathbf{M}\) 矩阵可能等于 \(diag(\mathbf{A})\) 或者容易求逆的 \(\mathbf{A}\) 矩阵的任意一种近似矩阵。

PCG 迭代的完整的伪代码如下 (这里原方程组设为 \(\mathbf{Ax} = \mathbf{b}\) ) $$ \[\begin{aligned} &\mathbf{x}=0 ; \mathbf{r}=\mathbf{b}-\mathbf{A x} ; \mathrm{k}=1 \\ & \text{do while(} \mathbf{r'r} \text{ not sufficiently small)} \\ & \quad \mathbf{z}=\mathbf{M}^{-1} \mathbf{r} \\ & \quad \tau_{\mathrm{k}-1}=\mathbf{z}^{\prime} \mathbf{r} \\ & \quad \text{if} (\mathrm{k}==1) \text{ then} \\ & \quad \quad \beta=0 ; \mathbf{p}=\mathbf{z} \\ & \quad \text{else} \\ & \quad \quad \beta=\tau_{\mathrm{k}-1} / \tau_{\mathrm{k}-2} ; \mathbf{p}=\mathbf{z}+\beta \mathbf{p} \\ & \quad \text { end if } \\ & \quad \mathbf{w}=\mathbf{A p} \\ & \quad \alpha=\tau_{\mathrm{k}-1} /\left(\mathbf{p}^{\prime} \mathbf{w}\right) \\ & \quad \mathbf{x}=\mathbf{x}+\alpha \mathbf{p} \\ & \quad \text { if }(\bmod (\mathrm{k}, 100) /=0) \text { then } \\ & \quad \quad \mathbf{r}=\mathbf{r}-\alpha \mathbf{w} \\ & \quad \text { else } \\ & \quad \quad \mathbf{r}=\mathbf{b}-\mathbf{A} \mathbf{x} \\ & \quad \text { end if } \\ & \quad \mathrm{k}=\mathrm{k}+1\\ &\text{end do} \\ \end{aligned}\]

$$ 尽管看上去比较复杂，但是 PCG 方法很容易应用，因为这里唯一的复杂计算就是 \(\mathbf{Ax}\) 和 \(\mathbf{Ap}\) 。

迭代的收敛性

“真实的”收敛标准为下式，其中 \(\mathbf{b}\) 为精确解。 \[ \frac{\|\mathbf{b^{(t)}-b}\|}{\|\mathbf{b}\|} \] 但是由于实际上我们不知道 \(\mathbf{b}\) ，因此在迭代过程中，我们常用下面的两个收敛标准：

第一种方式，根据两次迭代解的差距，公式如下。也就是说，如果两次迭代的解相差很小，则认为收敛。 \[ C_{d}^{(t)} = \frac{\|\mathbf{b^{(t)}-b^{(t-1)}}\|}{\|\mathbf{b^{(t)}}\|}=\frac{\sum_{i=1}^{p}\left(b_{i}^{(t)}-b_{i}^{(t-1)}\right)^{2}}{\sum_{i=1}^{p}\left(b_{i}^{(t)}\right)^{2}}<\varepsilon \] 第二种方式，根据左手项和右手项的差距，公式如下 \[ C_{r}^{(t)} = \frac{\|\mathbf{r-Ab^{(t)}}\|}{\|\mathbf{r}\|} < \varepsilon \] 其中 ，\(\varepsilon> 0\) 是预先给定的允许误差，通常是一个很小的数值。在每次迭代结束时，都检查上式是否成立，如果成立，就意味着本次迭代与上次迭代所得到的解已非常接近，可以认为迭代收敛，从而结束迭代，将当前迭代的解作为方程组的解。

混合模型方程组的迭代求解

直接迭代求解

直接迭代求解是先计算出方程组系数矩阵和右手项的各个元素，然后用前面所介绍的方法进行迭代求解，这种解法也称为对方程组迭代 (iteration on equations) 。一般来说，当数据可以读入内存中时，我们会用 GS 或 SOR 迭代法。但是，PCG 方法现在越来越流行，比如将预处理矩阵设为对角矩阵的 PCG 方法可以对所有的模型均收敛，包括复杂性状，随机回归模型，母性效应模型等；PCG 方法的收敛速度一般也优于前几种方法；PCG 方法还不需要确定任何参数。PCG 方法的一个劣势是它需要更多的内存。

间接迭代解法

用直接迭代求解并将方程组存储在硬盘上的缺陷在于从硬盘读取数据的次数太多，计算速度减慢，同时编程也较复杂。针对这种缺陷，Schaeffer 和 Kennedy (1986a, 1986b) 及 Misztal 和 Gianola (1987) 等分别提出了混合模型方程组的间接迭代解法，用这种方法不需要建立方程组，而是在每次迭代中读入观测值书数据文件和系谱数据文件，并同时计算方程组的解，故这种方法又称为对数据文件迭代 (iteration on data) 。因为在多数情况下混合模型方程组的大小要比数据文件的大小要大，这种方法所需的读取次数少于直接迭代方法，因而计算速度要快一点，这种方法的另一个优点是它的编程非常简单。