0先看目标：我们要读懂哪一行公式？

动物模型（Animal Model）几乎都写成下面这一行。本教程只讲读懂它所需的线性代数，不多也不少。

y = Xb + Zu + e，　　Var(u) = Aσa2，　Var(e) = Iσe2
y：所有个体的表型记录（一列数 → 向量）
b：固定效应，如性别、场、年季（向量）
u：每头动物的育种值（向量，我们最想要的）
e：残差（向量）
X、Z：把效应"分配"到记录上的 0/1 表（设计矩阵）
A：个体间的亲缘关系（分子血缘相关矩阵）

最终求解育种值靠的是 Henderson 的混合模型方程（MME）：

X′XX′Z
Z′XZ′Z + A−1α
b̂
û
 =
X′y
Z′y
 　α = σe2/σa2

看懂这一块需要：向量、矩阵、乘法、转置（′）、逆（−1）、方程组、对称/正定矩阵。下面一个个来。

贯穿全文的小例子

8 头犊牛，1–3 号是基础群（父母未知），4–8 号有断奶增重记录（kg）。这是 Mrode《Linear Models for the Prediction of Animal Breeding Values》第 3 章的经典例子，假设 σa2=20，σe2=40。

动物父母性别增重 y
1––––
2––––
3––––
41–公4.5
532母2.9
612母3.9
745公3.5
836公5.0
1 2 3 4 5 6 7 8 方框=公，圆=母（1、3 在系谱中作父亲）

1向量：一支箭头，也是一列数

向量就是"按顺序排好的一串数"。2 个数可以画成平面上的箭头，5 个数虽然画不出来，但规则完全一样。

a=(3,1) b=(1,2) a+b=(4,3)

加法：对应位置相加。几何上就是"箭头首尾相接"。
数乘：2a = (6, 2)，箭头方向不变、长度翻倍。

在动物模型里，y 就是一个 5 维列向量：

y =

点积（内积）：矩阵乘法的"零件"

两个等长向量对应相乘再相加：a·b = 3×1 + 1×2 = 5。几何含义：一个向量在另一个方向上的"投影"有多长。点积为 0 ⇔ 两箭头互相垂直。

记住这一点就够了

后面所有矩阵乘法，本质都是"行向量 · 列向量"的点积重复很多次。

2矩阵：一张表，也是一种"变形"

矩阵有两种看法：① 一张有行有列的数表（存数据）；② 一台把向量变成新向量的机器（做变换）。

大小用"行 × 列"描述：X 是 5×2（5 条记录、2 个性别），A 是 8×8（8 头动物）。

用 2×2 矩阵

ab
cd
变换平面时，第 1 列 (a,c) 是原来 x 轴单位箭头 î 的新位置，第 2 列 (b,d) 是 y 轴单位箭头 ĵ 的新位置。整个网格跟着它们一起变形。拖动滑块试试：

这个画布会在 第 6 节（行列式） 和 第 7 节（逆矩阵） 再用到：阴影方块的面积 = |行列式|；能"变回去"的变换才有逆矩阵。

3矩阵乘法：行 · 列

结果第 i 行第 j 列 = 左矩阵第 i 行 与 右矩阵第 j 列 的点积。把鼠标移到（或点击）右边结果矩阵的格子上看看它是怎么来的。

← 指向 C 中任一格

尺寸规则（最常出错）

(m×n) · (n×p) = m×p。
左边的"列数"必须等于右边的"行数"，中间两个数"抵消"。
例：X(5×2)·b(2×1) = 5×1，正好和 y 一样大。

顺序不能换

一般 AB ≠ BA，甚至 BA 可能尺寸不合根本算不了。几何上：先旋转再拉伸 ≠ 先拉伸再旋转。
但满足结合律：(AB)C = A(BC)。

矩阵乘向量 = 对矩阵的列做加权求和：Xb = b₁·(第 1 列) + b₂·(第 2 列)。这正是"每条记录加上它所属性别的效应"。

4转置与对称矩阵

转置：把行变成列（沿主对角线翻折）。育种文献常写 X′，数学书常写 XT，是同一个意思。

′ =
  • (AB)′ = B′A′（顺序反过来，像"穿衣脱衣"）
  • X′X 总是方阵，而且总是对称的

对称矩阵：A′ = A

以主对角线为镜子左右对称。动物模型里的 A、V、整个 MME 左边矩阵都是对称的——因为"牛 4 和牛 7 的亲缘"当然等于"牛 7 和牛 4 的亲缘"。软件只需存一半，省一半内存。

5特殊矩阵与设计矩阵 X、Z

单位矩阵 I


乘谁都不变：IA = A，相当于数字 1。几何：什么都不做。

对角矩阵


只有对角线非零。几何：沿坐标轴各自拉伸。X′X 在单因子时就是对角阵。

零矩阵 0、方阵

全是 0 的矩阵；行数=列数的叫方阵。只有方阵才谈得上行列式与逆。

设计（关联）矩阵：0/1 的"记账本"

每一行是一条记录，每一列是一个效应水平；该记录属于哪个水平，那一格就是 1。

X（5×2）：第 1 列"公"、第 2 列"母"。Z（5×8）：每列是一头动物，1–3 号没有记录所以整列为 0。

直觉：X′X 和 Z′Z 在"数个数"

X′X =　X′y =
X′X 对角线 = 每个性别有几条记录（公 3、母 2）；X′y = 每个性别的增重总和。所以 (X′X)−1X′y 就是"各性别的平均值"——这就是最小二乘。Z′Z 同理：对角线是每头动物的记录数（这里 1–3 号为 0，4–8 号为 1）。

6行列式：变换把面积放大了几倍

对 2×2 矩阵，det = ad − bc。几何意义：单位正方形变换后的（有方向的）面积。

  • |det| = 2：面积变 2 倍；det 为负：图形被"翻面"（镜像）。
  • det = 0：平面被压成一条线（或一个点）——信息丢失了，无法还原。回到第 2 节画布点"压扁（奇异）"看看。
  • 高维也一样：3×3 是体积缩放，n×n 是"n 维体积"缩放。

和育种的联系

det = 0 的矩阵叫奇异矩阵，没有逆。MME 左边矩阵如果奇异（比如模型里同时放了总均值和全部性别水平），就需要第 8 节的处理。REML 估计方差组分时还会用到 log|V| 或 log|C|。

7逆矩阵与线性方程组

逆矩阵 A−1 就是"撤销"变换：A−1A = I。它是矩阵世界的"除法"。

2×2 求逆公式

ab
cd
−1 = 1/(ad−bc) ×
d−b
−ca
主对角线互换、副对角线变号，再除以行列式——这就是为什么 det=0 时不能求逆（除以 0）。

方程组 = 找交点

方程组 Ax = b 的解是 x = A−1b。两个未知数时，每个方程是一条直线，解就是交点。固定第一个方程 x + y = 4，调节第二个方程：

第二个方程：a·x + b·y = c

MME 就是这样一个方程组

只不过未知数有 2 + 8 = 10 个（2 个性别效应 + 8 个育种值），"直线"变成了 10 维空间里的"超平面"。原理完全相同：系数矩阵行列式不为 0 → 唯一解。大型评估（几百万头）不会真去求逆，而是用迭代法（如 Gauss-Seidel、PCG）直接解方程。

8秩、线性相关与广义逆

秩（rank）= 矩阵里真正"独立"的列（或行）有几个 = 变换后空间还剩几维。

如果模型里同时放总均值 μ 和两个性别水平：

X =　X′X =

第 1 列 = 第 2 列 + 第 3 列，有一列是"多余的"（线性相关），秩 = 2 而不是 3，所以 det(X′X) = 0。直观上：给 μ 加 10、给两个性别各减 10，预测值完全不变——数据无法区分，解有无穷多个。

三种常见处理

① 去掉一列（如 SAS 把最后一个水平设为 0）；② 加约束（如两性别效应之和为 0）；③ 用广义逆 G（满足 AGA = A）得到一组解。
关键结论：虽然 b̂ 本身不唯一，但可估函数（如"公 − 母"的差）以及育种值 û 的排序是唯一的。所以育种值通常以"相对值"解读。

9方差协方差矩阵：一个"椭圆"

两个随机变量的方差和协方差可以排成对称矩阵 V = [[Var₁, Cov],[Cov, Var₂]]。数据点云的形状就是一个椭圆。

  • Cov = 0：椭圆的轴与坐标轴平行；Cov > 0：向右上倾斜（正相关）。
  • 椭圆的长短轴方向 = 特征向量；轴长² ∝ 特征值。
  • 把 Cov 拉得太大（|Cov| > √(Var₁·Var₂)，即相关系数 |r|>1），就会出现负特征值——这样的"方差矩阵"在现实中不可能存在。

正定矩阵（positive definite）

所有特征值 > 0，等价于对任意非零 x 都有 x′Vx > 0（"任何组合的方差都为正"）。合法的方差协方差矩阵必须正定（或至少半正定）。A、G、R 都应如此；多性状分析中估出的遗传协方差矩阵"不正定"是常见报错，原因往往就是相关系数估出了界。

多性状模型里的克罗内克积 ⊗

两个性状的育种值方差：Var(u) = G₀ ⊗ A。⊗ 的规则：用左边矩阵的每个元素去乘整个右边矩阵，拼成大块。

G₀ = 2 个性状的遗传方差/协方差，A = 2 头半同胞（亲缘 0.5）。结果 4×4 描述"任意两头动物 × 任意两个性状"的协方差。

10亲缘关系矩阵 A（分子血缘相关矩阵）

aij = 动物 i 与 j 的加性遗传相关的 2 倍亲缘系数；aii = 1 + Fi（F 为近交系数）。

表格法（按出生顺序，父母在前）

设动物 j 的父母为 s、d（未知则视为 0）：
・非对角：aij = ½ (ais + aid)　——"我和你的亲缘 = 我和你父母亲缘的平均"
・对角：ajj = 1 + ½ asd　——"父母越亲，我越近交"
A =

读读看：a14 = 0.5（父子）；a56 = 0.25（同母半同胞）；8 号的父母 3 和 6 之间 a36 = 0，所以 8 号不近交，a88 = 1。

为什么 MME 用 A−1 而不是 A

A 几乎处处非零（稠密），而 A−1 非常稀疏：只有"父母–子女"和"父–母"之间才非零，并且可以直接从系谱写出来（Henderson 规则），根本不用先算 A 再求逆。这就是几百万头动物也能做 BLUP 的秘密。
A−1 =

11从模型推导方差：V、协方差与 BLUP 公式

模型 y = Xb + Zu + e 只写了"均值结构"。把几条假设和三条运算规则用上，就能推出数据的方差结构，进而推出育种值该怎么估。

第 1 步：模型假设（推导的全部"原料"）

假设含义生物学解释
E(u) = 0，E(e) = 0随机效应均值为 0育种值是相对于群体的离差；"平均"部分已放进 Xb
Var(u) = Aσₐ² ＝ G育种值的协方差由亲缘决定只考虑加性效应；亲缘越近，育种值越相像
Var(e) = Iσₑ² ＝ R残差互不相关、方差相同没有共同环境（如同窝效应）；若有，需另加随机效应
Cov(u, e′) = 0遗传与环境无关好基因的牛不会被系统地分到好圈舍
b 固定、σₐ²、σₑ² 已知BLUP 需要先知道方差组分实际中先用 REML 估计（见本节末）

第 2 步：三条运算规则（标量 → 矩阵）

① 期望：线性的

E(Ax + c) = AE(x) + c
标量版：E(ax+c) = aE(x)+c

② 方差："夹心"

Var(Ax) = A·Var(x)·A′
标量版：Var(ax) = a²Var(x)。a² 拆成"左乘 A、右乘 A′"。几何上就是：把第 9 节的椭圆用第 2 节的矩阵 A 变形一下。

③ 和的方差

Var(x+w) = Var(x) + Var(w) + Cov(x,w′) + Cov(w,x′)
Cov(Ax, (Bw)′) = ACov(x,w′)B′

第 3 步：推导

均值：E(y) = Xb + ZE(u) + E(e) = Xb

方差：Xb 是常数，不贡献方差，所以

Var(y) = Var(Zu + e) 　(规则③，Cov(u,e′)=0 两项消失)
= Z·Var(u)·Z′ + Var(e) 　(规则②)
= ZAZ′σₐ² + Iσₑ² ＝ V

记录与育种值的协方差：

Cov(y, u′) = Cov(Zu + e, u′) = Z·Var(u) = ZAσₐ²

把例子代进去（σₐ²=20，σₑ²=40），V 的行列对应有记录的动物 4–8：

V =（行列 = 动物 4,5,6,7,8）

看懂 V 的每个元素，就是经典数量遗传学

・对角线 Vii = aiiσₐ² + σₑ² = 20 + 40 = 60 = 表型方差 σₚ²，h² = 20/60 = 1/3。
・非对角 Vij = aijσₐ²：
　这正是教科书上"半同胞表型协方差 = ¼σₐ²、亲子协方差 = ½σₐ²"的矩阵写法。
・亲属之间表型相关，所以 y 不能当作独立数据处理——这就是普通最小二乘不够用的原因。

第 4 步：从方差到"最好的预测"

已知 y 后，怎样猜 u 最准？统计上的答案是条件期望。在正态假设下，它是一条回归：

û = Cov(u, y′) · Var(y)−1 · (y − E(y)) = GZ′V−1(y − Xb̂)

其中 b̂ 用广义最小二乘（GLS）估计：b̂ = (X′V−1X)−1X′V−1y。V−1 就像"权重"：方差大、彼此相关的记录被打折。

标量直觉：只有自身一条记录时

所有矩阵都变成数：G = σₐ²，Z = 1，V = σₚ²，于是
û = (σₐ² / σₚ²)(y − μ) = h²(y − μ)
这就是经典的"育种值 = 遗传力 × 表型离差"。回归系数 = 协方差 / 方差。矩阵版 GZ′V−1 只是把它推广到"多条记录、亲属互相借信息"的情况。拖动滑块看：点是模拟的个体，横轴表型离差，纵轴真实育种值。
红色虚线：斜率 1，直接拿表型当育种值。
绿色实线：斜率 h²，即 BLUP 的"收缩"预测，离差平方和最小。
h² 越低，点云越扁平，最佳预测线越平——证据少，就越保守。

Henderson 的贡献：不用求 V−1

V 是 n×n 的稠密矩阵（n = 记录数，可达上亿），求逆不现实。Henderson 证明：上面 GLS + 回归公式的解，与混合模型方程的解完全相同；而 MME 只用到稀疏的 A−1 和 R−1。下面在浏览器里用两种方法分别算一遍：

第 5 步：预测有多准——PEV 与可靠性

预测误差 u − û 的方差（PEV）可以直接从 MME 左边矩阵的逆 C−1 读出：

Var(ui − ûi) = ciiσₑ²，　　r²i = 1 − PEVi / (aiiσₐ²)，　　SEP = √PEV

其中 cii 是 C−1 里动物 i 对应的对角元。还有两个常用结论：Var(û) = Cov(u, û′) = G − PEV（估计值的方差比真实育种值小，也就是"收缩"），以及 r² = Var(ûi)/Var(ui)。

附：σₐ²、σₑ² 本身怎么来？（REML 一瞥）

以上都假设方差组分已知。实际中先用 REML（约束最大似然）估计：它去掉固定效应 b 的影响，再让数据出现的"可能性"最大。其对数似然为
log L = −½ [ log|V| + log|X′V−1X| + y′Py ] + 常数，　P = V−1 − V−1X(X′V−1X)−1X′V−1
这里用到了行列式（第 6 节）、逆（第 7 节）。EM-REML 的迭代式更直观（q 为动物数，Cuu 为 C−1 中育种值那一块）：
σ̂ₐ² = [ û′A−1û + tr(A−1Cuu)σₑ² ] / q
直观上：第一项是"估出来的育种值有多分散"，第二项把预测误差补回来，因为 û 被收缩了。这里又出现了迹 tr（对角线之和）。

12把一切拼起来：混合模型方程（MME）

下面所有数字都由浏览器实时计算。拖动遗传力，看育种值如何变化。

第 1 步：拼出左边矩阵 C 与右边向量

蓝 = X′X，绿 = X′Z 和 Z′X，红 = Z′Z + A−1α。注意它是对称的。

C =　r =

第 2 步：解方程 C·[b̂; û] = r

[b̂; û] = C−1r =
估计育种值 EBV（kg）

你能从结果里读出什么

・1–3 号没有记录，但通过 A−1 从后代那里"借"到了信息，所以也有育种值。
・α 是"收缩"力度：h² 越低，α 越大，A−1α 在对角线上越重，育种值越被拉向 0（越保守）。把 h² 拖到 0.05 看看，再拖到 0.9。
・Z′Z + A−1α 中的"+"让矩阵变得可逆：即使有动物没记录（Z′Z 对角为 0），加上 A−1α 后仍正定。
・预测误差方差 PEV 和可靠性也来自 C−1（推导见第 11 节）：PEV = ciiσₑ²，r² = 1 − PEV/(aiiσₐ²)；本例无近交（aii=1），化简为 r² = 1 − cii·α。

13速查表

符号名称一句话在动物模型里
a·b点积对应相乘再相加矩阵乘法的每一格
AB矩阵乘法行·列；(m×n)(n×p)=m×pXb、Zu、X′X …
A′ / Aᵀ转置行列互换；(AB)′=B′A′X′y、Z′Z
I单位矩阵相当于 1，什么都不做Var(e)=Iσₑ²
|A|, det行列式面积/体积缩放倍数，0 ⇒ 不可逆REML 似然中的 log|V|
A⁻¹逆矩阵撤销变换；A⁻¹A=IA⁻¹α 加进 MME
rank秩独立列的个数固定效应不满秩 ⇒ 广义逆
A⁻广义逆AGA=A，奇异时也能"凑"一个解b̂ 不唯一，可估函数唯一
λ, v特征值/向量Av=λv，椭圆的轴判断正定、主成分
正定positive definite特征值全 >0，x′Vx>0合法的 A、G、R
⊗克罗内克积每个元素乘整个矩阵多性状 G₀⊗A
tr(A)迹对角线之和REML/EM 算法公式

进一步阅读：Mrode R.A.《Linear Models for the Prediction of Animal Breeding Values》第 1–3 章；张沅《家畜育种学》中 BLUP 部分。有了本页的基础，读这些书时遇到的矩阵运算应当都不陌生了。