上周帮 A 同学看一个回归模型，他把“渠道类型”这一列直接按 1、2、3 塞进去跑线性回归，跑完问我：线上渠道的系数是 -0.4，是不是说明线上效果差？我问他 1、2、3 分别代表什么，他说 1 是自然流量、2 是投放、3 是老带新。问题到这儿就清楚了——他把一个分类标签当成了连续强度，模型算出来的斜率没有任何业务含义，只是“标签数字每加 1”这个动作的数学后果。这类事我在项目里见得太多了，几乎每隔一段时间就要重新解释一遍：哑变量到底在做什么、为什么要留一个参照组、为什么全部放进去软件就报共线、为什么换了参照组之后所有系数都变了但预测值一个没变。这一串问题背后其实是同一个东西——参数可识别性。哑变量不只是一个编码技巧，它是一套关于“如何在不可逆的标签世界里建立坐标系”的数学约定；参照系选择决定了参数的语义，而可识别性决定了你能不能从数据里唯一地把这些参数解出来。这篇就把这条线索从头到尾捋一遍，公式、数字、代码、坑，都摆出来。适合已经会用回归但被系数解释折磨过的人，也适合刚接触特征工程、想知道get_dummies背后到底发生了什么的同学。

1. 哑变量在数学上到底做了什么

1.1 把标签翻译成坐标：一个直观的类比

分类变量天生不具备数序结构。“自然流量”“投放”“老带新”这三个词之间没有大小、没有距离、没有倍数关系，你硬给它赋 1、2、3，等于强行宣布“投放是自然流量的两倍”“老带新是自然流量的三倍”，这个假设没有任何业务依据支撑，模型却会老老实实按这个假设去拟合，最后系数自然没法解释。

哑变量的思路是换一套表达：不给每个类别一个“数值”，而是给每个类别一个“位置”。就像在一个三维空间里，三个类别分别占据三个坐标轴方向上的单位向量——A 是 (1,0,0)，B 是 (0,1,0)，C 是 (0,0,1)。这样任意两个类别之间的距离都相同，任意两个类别之间都不存在“谁比谁大”的暗示，模型只需要为每个方向学一个偏移量，就能拟合出各组的不同水平。这是哑变量最本质的东西：它把“分类”这个定性概念，翻译成了“互斥指示”这个定量结构。

但这里立刻出现一个麻烦。如果三个类别各占一个方向，再加上一个全局截距（可以理解为原点位置），那这个坐标系就多了一个自由度——原点和三个方向之间存在一组线性关系：三个方向向量加起来正好等于全 1 向量，而全 1 向量就是截距列。换句话说，你既指定了原点，又指定了三个方向，但其中一个方向可以由另外两个方向加原点推出来，信息是重复的。数学上这叫共线，工程上这叫冗余，统计上这叫不可识别。参照系的选择，本质上就是在处理这个冗余——你必须主动放弃一个方向，让剩下的方向相对于被放弃的那个方向来定位。

1.2 形式化定义与设计矩阵

把上面的话写成公式。设分类变量 $X$ 取值于 ${1,2,\dots,K}$，共 $K$ 个水平。定义一组指示函数：

$$D_k^{(i)} = \begin{cases} 1, & X^{(i)} = k \ 0, & X^{(i)} \neq k \end{cases}$$

对第 $i$ 个样本，$K$ 个指示变量里恰好有一个等于 1，其余全为 0。这就是哑变量的定义，也叫 one-hot 编码。注意它是一个互斥且完备的编码：互斥意味着不会有两个同时为 1，完备意味着不会全为 0（前提是这 $K$ 个水平覆盖了所有可能取值，没有遗漏、没有空组）。

把样本堆起来，就得到设计矩阵。以 $K=3$、样本量 $n$ 为例，完整写法是：

$$X_{\text{full}} = \begin{bmatrix} 1 & D_1^{(1)} & D_2^{(1)} & D_3^{(1)} \ 1 & D_1^{(2)} & D_2^{(2)} & D_3^{(2)} \ \vdots & \vdots & \vdots & \vdots \ 1 & D_1^{(n)} & D_2^{(n)} & D_3^{(n)} \end{bmatrix}$$

这个矩阵有 $n$ 行 4 列。第 1 列是全 1 的截距列，后 3 列是哑变量。行方向看，每一行形如 $(1,1,0,0)$ 或 $(1,0,1,0)$ 或 $(1,0,0,1)$，取决于该样本属于哪一组。注意每一行的后三列之和恒等于 1，也就是：

$$D_1 + D_2 + D_3 = \mathbf{1}$$

这个恒等式看起来平平无奇，但它是后面所有麻烦的总源头。它意味着后三列的和等于第一列，也就是说这四列线性相关，矩阵的秩最多是 3，而列数是 4——秩亏 1。用线性代数的语言，$X_{\text{full}}$ 的零空间至少包含向量 $(0,1,1,1)^\top$ 以及它的任意倍数，因为：

$$X_{\text{full}} \cdot (0,1,1,1)^\top = 0 \cdot \mathbf{1} + D_1 + D_2 + D_3 = \mathbf{1} \neq 0$$

等一下，这里要小心，符号写反了。正确的零向量是 $(1,-1,-1,-1)^\top$：

$$X_{\text{full}} \cdot (1,-1,-1,-1)^\top = 1\cdot\mathbf{1} - D_1 - D_2 - D_3 = \mathbf{1} - \mathbf{1} = \mathbf{0}$$

对，这样才等于零向量。这个 $(1,-1,-1,-1)^\top$ 就是那个“看不见的冗余方向”：如果你把截距加上 1、同时把三个哑变量的系数各减 1，模型对所有样本的预测值完全不变。参数变了，预测没变——这就是不可识别的典型症状。

1.3 参照系是哪一列，为什么它“消失”了

工程上最常见的处理是删掉一列。删掉 $D_3$ 之后，设计矩阵变成 $n \times 3$：

$$X = [\mathbf{1}, D_1, D_2]$$

现在三列线性无关（前提是每组至少有一个样本，这点后面会说），秩为 3，列满秩，参数可以唯一解出。被删掉的 $D_3$ 对应的类别就叫参照组（reference level），也常被叫作基线组、对照组。它并没有真的“消失”——它的信息被编码进了截距里。

这一点很多人第一次接触时会绕不过来。做个具体的推导就明白了。假设模型里只有这一个分类变量，那么第 $k$ 组的预测均值是：

• 参照组（$k=3$）：$\hat{y} = \beta_0$
• 第 1 组：$\hat{y} = \beta_0 + \beta_1$
• 第 2 组：$\hat{y} = \beta_0 + \beta_2$

用最小二乘解出来，$\hat\beta_0$ 就是第 3 组的样本均值，$\hat\beta_1$ 是第 1 组均值减第 3 组均值，$\hat\beta_2$ 是第 2 组均值减第 3 组均值。所以：

被删掉的那一列不是不存在，而是它“化成了原点”。截距吸收了参照组的水平，其余系数全部变成相对于参照组的差值。

这个理解方式很有用，因为它一次性解释了三件事：为什么截距不能随便解释成“所有组的平均值”（它只是参照组的平均值）；为什么哑变量系数是“差值”而不是“水平”；为什么换一个参照组，截距和所有哑变量系数都会变，但组与组之间的差值不会变。

顺带说一句，有些教材会把参照组说成“基准类别”，容易让人误以为它更重要或者更特殊。它一点都不特殊，纯粹是数学上的锚点。选谁当锚点不改变模型对数据的拟合能力，只改变你读系数时的参照对象。

2. 参照系选择的连锁反应：截距、系数与解释口径

2.1 三种等价参数化与系数换算公式

同一个分类变量，至少有三种常见的满秩参数化方式，它们生成的列空间完全相同，因此拟合值、残差、残差平方和、$R^2$（在有截距定义下）全部一样，但系数含义截然不同。

第一种：截距 + $K-1$ 个哑变量。最常见，也是 R 的lm、Python 的statsmodels默认行为（把第一个水平当参照）。系数解释为“该组与参照组的均值差”。

第二种：无截距 + $K$ 个哑变量。去掉全 1 列，保留所有哑变量。此时每一列直接对应一组的均值，系数就是该组均值本身，读起来非常直观。但它有两处副作用需要留意：一是模型的 $R^2$ 定义会变（后面细讲），二是如果某个水平在数据里没有样本，那一列全 0，矩阵立刻秩亏。

第三种：效应编码（effect coding / sum-to-zero）。保留截距，但把参照组那一列的编码从 0 改成 -1。也就是对非参照组，$D_k$ 取 1 或 0；对参照组，全部 $K-1$ 个哑变量都取 -1。这样每一行的哑变量之和不再恒等于某个常数，共线关系被打断。系数解释为“该组与各组均值（简单平均）的偏差”。

第三种最容易被误解，因为它的截距不是加权总均值。举个例子：三组样本量分别是 3、2、5，组均值分别是 10、12、15。加权总均值是 $(3\times10 + 2\times12 + 5\times15)/10 = 12.9$，但效应编码解出来的截距是 $(10+12+15)/3 = 12.333$。它对应的是“各组均值再取平均”，也就是假设三组样本量相等时的总均值。这个细节在样本不均衡的时候会造成解释偏差，写报告时如果直接说“总体平均水平是 12.33”，而实际样本均值是 12.9，就容易被质疑。用之前先确认你要的到底是哪一个。

三种参数化之间的换算有固定公式。设原来以第 $K$ 组为参照，截距为 $\beta_0$，系数为 $\beta_1,\dots,\beta_{K-1}$。现在想把参照组换成第 $j$ 组（$j < K$），新参数记作 $\beta_0', \beta_k'$。由于两套参数必须给出相同的组均值，逐组列方程就能解出来：

$$\beta_0' = \beta_0 + \beta_j, \qquad \beta_k' = \beta_k - \beta_j \ (k \neq j),\qquad \beta_K' = -\beta_j$$

验证一下逻辑：原来第 $j$ 组的均值是 $\beta_0+\beta_j$，新模型里 $j$ 成了参照组，均值应该等于新截距 $\beta_0'$，所以 $\beta_0' = \beta_0+\beta_j$，成立。原来第 $K$ 组的均值是 $\beta_0$，新模型里它变成了普通组，均值是 $\beta_0' + \beta_K'$，代入得 $\beta_0+\beta_j+\beta_K' = \beta_0$，所以 $\beta_K' = -\beta_j$，也成立。

这个公式不只是用来做手工换算，它还解释了一个常见现象：换了参照组之后，某些系数会“变号”。因为 $\beta_K' = -\beta_j$，如果原来第 $j$ 组显著高于参照组（$\beta_j > 0$），那么新模型里原来的参照组相对新参照组就是负的。这完全是坐标系翻转的结果，不是数据变了，也不是模型变差了。

2.2 用一组真实数字走一遍

光看公式容易飘，拿一组具体数字落地。还是上面那三组数据：

总样本量 $n=10$，总和 129，加权均值 12.9。组内残差平方和：A 组 $(9-10)^2+(10-10)^2+(11-10)^2=2$，B 组和 C 组都是 0，所以 $SSE=2$。总平方和 $SST = 50.9$，$R^2 = 1 - 2/50.9 \approx 0.9607$。

现在看四种参数化下的系数：

四种哑变量写法给出的拟合值完全一致，$SSE$ 都是 2。参照 C 时，A 的系数 -5 意思是“A 组比 C 组低 5 个单位”；参照 A 时，C 的系数 5 意思是“C 组比 A 组高 5 个单位”。同一个事实，两种表述，符号相反，都对。

效应编码那一行的截距 12.333 是三个组均值的简单平均，A 的 -2.333 表示 A 组比这个“中心”低 2.333。三组偏差之和为 $(-2.333) + (-0.333) + 2.667 = 0$，这是效应编码的固有性质，可以用来自查算得对不对。

再看一个容易被忽视的差异：无截距模型的 $R^2$ 不能和上面几个直接比。有截距时 $R^2 = 1 - SSE/SST$，其中 $SST$ 是围绕均值的离差平方和，得 0.9607。无截距时很多软件报的是未中心化的 $R^2 = 1 - SSE/\sum y_i^2$，算出来是 $1 - 2/1715 \approx 0.9988$。看着漂亮多了，但这个数字不能拿来说“模型更好了”，因为分母换了，衡量基准不一样。我见过有人为了“提高 $R^2$”把截距去掉，这纯粹是自欺欺人。

2.3 参照组怎么选才不给自己挖坑

参照组的选择不改变模型，但改变你向别人解释模型时的难度。我一般按三个优先级来定。

第一优先级：业务上的基线。如果存在一个天然的对照（未干预组、默认配置、老版本、空白对照），就选它。汇报时“相比默认配置，新策略提升了 X”这种句式，听众一秒就懂。反过来，如果软件按字母序默认选了“A 组”当参照，而 A 组其实是个只有 3 个样本的试验组，那你后面所有系数都是“相对于一个边缘组”的差值，解释起来非常别扭。

第二优先级：样本量足够大的组。这不只是心理安慰，有实打实的方差依据。在单因子模型里，可以精确推出：

$$\mathrm{Var}(\hat\beta_0) = \frac{\sigma^2}{n_{\text{ref}}}, \qquad \mathrm{Var}(\hat\beta_k) = \sigma^2\left(\frac{1}{n_k} + \frac{1}{n_{\text{ref}}}\right)$$

推导不复杂：设计矩阵 $X=[\mathbf{1}, D_1,\dots,D_{K-1}]$，其 Gram 矩阵分块为 $\begin{bmatrix} n & u^\top \ u & B\end{bmatrix}$，其中 $B=\mathrm{diag}(n_1,\dots,n_{K-1})$，$u=(n_1,\dots,n_{K-1})^\top$。用分块求逆公式算 $(1,1)$ 位置，得到 $1/(n - u^\top B^{-1}u) = 1/(n - \sum_{k<K} n_k) = 1/n_K$，也就是参照组样本量的倒数。有了这个，$\mathrm{Var}(\hat\beta_k)$ 就是把“第 $k$ 组均值减参照组均值”的方差直接相加，得到上面的式子。

代入刚才的数据，$\hat\sigma^2 = 2/7 \approx 0.2857$。参照 C（$n=5$）时，$\mathrm{SE}(\hat\beta_0)=\sqrt{0.2857/5}\approx 0.239$；参照 B（$n=2$）时，$\mathrm{SE}(\hat\beta_0)=\sqrt{0.2857/2}\approx 0.378$。同一份数据，截距的标准误差了将近 60%。

这里必须补一句严谨的话，否则容易过度推广：任意两组之间的直接对比，其方差 $\sigma^2(1/n_i+1/n_j)$ 与参照组选谁无关。参照组只影响“以参照组为锚点”的那些对比的方差。所以正确表述是：选大组做参照，能让截距和所有“相对参照”的系数更稳；但如果你真正关心的是 A 和 B 的差异，选谁当参照都不会让这个差异估得更准，该做的是直接构造对比或者用效应编码。

第三优先级：样本量为零的组绝对不能当参照。这条是硬约束。如果参照组一个样本都没有，$1/n_{\text{ref}}$ 趋于无穷，参数根本估不出来，软件会报错或者给出极端值。空组问题后面还会展开。

3. 参数可识别性：秩亏的数学本质

3.1 可识别性的定义：从“两个参数给出同一个模型”说起

可识别性（identifiability）听起来很学术，但它的定义非常朴素：如果两个不同的参数值能产生完全相同的观测分布，那这组参数就不可识别。注意这里的“不同参数值”和“相同分布”是针对整个模型族说的，和数据量无关——哪怕你有一亿条数据，只要模型本身存在这种等价关系，你也永远分不清到底是哪个参数值在起作用。

在线性模型 $E[y \mid X] = X\theta$ 里，这个定义可以化成一个纯代数条件。假设 $\theta_1 \neq \theta_2$ 但 $X\theta_1 = X\theta_2$，那么 $X(\theta_1-\theta_2)=0$，即 $\theta_1-\theta_2$ 落在 $X$ 的零空间里。所以：

$\theta$ 可识别 $\iff$ $X$ 的零空间只有零向量 $\iff$ $X$ 列满秩。

这个等价关系是整个话题的枢纽。它把“统计上能不能分清”翻译成了“矩阵的秩够不够”，而秩是可以直接算的。反过来，只要 $X$ 的零空间非平凡，就存在无穷多组参数给出完全相同的拟合结果——这就是不可识别。

回到虚拟变量陷阱。$X_{\text{full}} = [\mathbf{1}, D_1, D_2, D_3]$ 有 4 列，但 $D_1+D_2+D_3=\mathbf{1}$，所以 $\mathrm{rank}(X_{\text{full}}) \le 3$，零空间至少是一维的，由 $(1,-1,-1,-1)^\top$ 张成。这意味着：对任意实数 $t$，参数向量

$$\theta + t\cdot(1,-1,-1,-1)^\top = (\beta_0+t,\ \beta_1-t,\ \beta_2-t,\ \beta_3-t)^\top$$

给出的预测值完全相同，因为 $X_{\text{full}}(\theta + tv) = X_{\text{full}}\theta + tX_{\text{full}}v = X_{\text{full}}\theta$。$t$ 取什么值都行，无穷多个解。这就是“不可识别”的完整数学图像：不是解不出来，而是解太多。

顺便说，这个视角也解释了为什么“多加几个变量总没坏处”是一种危险的直觉。每加一个与其他列线性相关的变量，零空间就多一个维度，不可识别的方向就多一个。

3.2 虚拟变量陷阱的秩证明

把上面的论证写严谨一点，因为这是最常被问到的地方。

设 $X_{\text{full}}$ 是 $n\times (K+1)$ 矩阵，第一列为 $\mathbf{1}$，后面 $K$ 列为 $D_1,\dots,D_K$。对每一行，$D_1^{(i)}+\cdots+D_K^{(i)}=1$，所以按列求和：

$$D_1 + D_2 + \cdots + D_K = \mathbf{1}$$

这说明第 1 列（截距）可以写成后面 $K$ 列的线性组合，系数全为 1。因此这 $K+1$ 个列向量线性相关，$\mathrm{rank}(X_{\text{full}}) \le K$。

另一方面，如果每个水平至少有一个样本，那么 $D_1,\dots,D_K$ 这 $K$ 列本身是线性无关的（假设 $\sum c_k D_k = 0$，取第 $i$ 个样本，它只属于某一组 $k_0$，则 $c_{k_0}=0$；对每组各取一个样本，就得到所有 $c_k=0$），所以 $\mathrm{rank}(X_{\text{full}}) = K$，秩亏恰好是 1。

如果还有某个水平样本量为零，比如 $n_3 = 0$，那 $D_3$ 这一列全为 0，秩再降 1，变成 $K-1$，秩亏 2。这就是为什么“空组”比“虚拟变量陷阱”更隐蔽——它不会有那个漂亮的 $\sum D_k = \mathbf{1}$ 关系来提醒你，只是安静地把秩悄悄拿掉一格。

删掉任意一个哑变量之后，比如删掉 $D_K$，剩下 $[\mathbf{1}, D_1,\dots,D_{K-1}]$ 共 $K$ 列。这时 $\mathbf{1}$ 可以写成 $D_1+\cdots+D_K$，但 $D_K$ 不在列集合里了，所以 $\mathbf{1}$ 不再能被剩下的列线性表出，秩恢复到 $K$，列满秩，可识别。这就是“$K$ 个类别用 $K-1$ 个哑变量”这条经验的严格来源。

有一个细节值得强调：“秩亏”和“软件报错”不是一回事。R 的lm遇到共线会给出NA系数并在警告里说明有几个系数因奇异性未被定义，这是软件替你做了参数选择；而 NumPy 的lstsq、scikit-learn 的很多估计器会默默返回一个最小范数解，一声不吭地把数算出来。后者更危险，因为你拿到的系数看着正常，实际只是无穷多解里随便挑的一个，解释全靠运气。我的习惯是建模前先跑一遍秩检查，不依赖软件的报错。

3.3 可识别与可估是两件事

这里要引入一个更精细的概念：可估函数（estimable function）。参数不可识别，不代表所有关于参数的问题都回答不了。

定义是：线性组合 $c^\top\theta$ 可估，当且仅当 $c$ 落在 $X^\top X$ 的行空间里（等价地，$c$ 与 $X$ 的零空间正交）。直观理解：如果两个参数值给出相同的预测，那么任何“不区分这两个参数值”的组合就是可估的。

在虚拟变量陷阱那个例子里，$\theta = (\beta_0,\beta_1,\beta_2,\beta_3)$ 不可识别，零空间方向是 $v=(1,-1,-1,-1)^\top$。那么一个线性组合 $c^\top\theta$ 可估，需要 $c^\top v = 0$，即 $c_0 - c_1 - c_2 - c_3 = 0$。检查几个有意思的组合：

• 组 A 的均值 $\beta_0+\beta_1$，对应 $c=(1,1,0,0)$，$c^\top v = 1-1 = 0$，可估。
• 组 A 与组 B 的差 $\beta_1-\beta_2$，对应 $c=(0,1,-1,0)$，$c^\top v = 0-1+1=0$，可估。
• 组 A 与参照组 C 的差 $\beta_1-\beta_3$，对应 $c=(0,1,0,-1)$，$c^\top v = -1+1 = 0$，可估。
• 单个 $\beta_0$ 本身，对应 $c=(1,0,0,0)$，$c^\top v = 1 \neq 0$，不可估。

结论很清晰：即使你误放了全部哑变量加截距，组均值、组间差这些你真正关心的量依然可以被唯一估计出来，只是单个参数没有意义。这也解释了为什么某些软件在秩亏时仍然能给出正确的预测和正确的组间对比检验——它内部用的是广义逆或者约束优化，绕开了不可识别的部分。

但这里有个陷阱。不同的广义逆给出的特解不同。$X^\top X$ 的 Moore-Penrose 逆给出的最小范数解，恰好是“让参数向量长度最小的那个解”，这个选择没有任何业务含义。所以如果你用lstsq拿到一组数就去解释，很可能解释的是一个数学约定，而不是数据事实。判断方法很简单：换一种编码方式重跑一遍，如果某个系数变了，那它就不是可估的，别解释。

4. 实操：构造矩阵、诊断秩亏、验证等价性

4.1 手写三种编码并比对秩与条件数

理论说完，动手跑一遍。下面这段代码用 NumPy 和 pandas 手搓三种编码，打印秩和条件数。代码不长，但每一步都对应上面讲的一个概念。

import numpy as np import pandas as pd # 三组，样本量不等 y = np.array([9., 10., 11., 12., 12., 15., 15., 15., 15., 15.]) g = np.array(['A']*3 + ['B']*2 + ['C']*5) n = len(y) D = pd.get_dummies(pd.Series(g), prefix='g').astype(float) print('水平:', D.columns.tolist()) # ['g_A', 'g_B', 'g_C'] # 方案一：截距 + K-1 个哑变量（参照 C） X1 = np.column_stack([np.ones(n), D[['g_A', 'g_B']].values]) # 方案二：无截距 + K 个哑变量 X2 = D.values # 方案三：截距 + 全部 K 个哑变量（陷阱） X3 = np.column_stack([np.ones(n), D.values]) def report(name, X): s = np.linalg.svd(X, compute_uv=False) r = np.linalg.matrix_rank(X) cond = s[0] / s[-1] if s[-1] > 1e-10 else float('inf') print(f'{name:16s} 列数={X.shape[1]} 秩={r} 最小奇异值={s[-1]:.2e} 条件数={cond:.3e}') report('截距+K-1', X1) report('无截距+K', X2) report('截距+K(陷阱)', X3)

跑出来的结果大致是：方案一列数 3、秩 3、条件数在个位数；方案二列数 3、秩 3、条件数也不大；方案三列数 4、秩 3、最小奇异值约 $10^{-16}$ 量级、条件数报inf。

这里要特别注意那个最小奇异值。理论上它严格等于 0，但浮点运算下会是一个极小的正数，比如 $6.7\times10^{-16}$。如果你用s[-1] > 0来判断是否共线，或者用s[0]/s[-1]算条件数而不设阈值，会得到一个巨大的有限值（$10^{16}$ 量级）而不是inf。工程上的做法是设一个相对阈值：s[-1] < 1e-10 * s[0]就判定为秩亏。这个阈值不是随便定的，$10^{-10}$ 大致相当于双精度浮点数在 $n$ 规模矩阵上的舍入误差累积上限，比它小的奇异值基本可以认为是数值噪声。

再看看系数，验证“列空间相同则拟合值相同”：

b1, *_ = np.linalg.lstsq(X1, y, rcond=None) b2, *_ = np.linalg.lstsq(X2, y, rcond=None) b3, *_ = np.linalg.lstsq(X3, y, rcond=None) print('参照C :', np.round(b1, 3)) # [15. -5. -3.] print('无截距:', np.round(b2, 3)) # [10. 12. 15.] print('陷阱 :', np.round(b3, 3)) # 某个最小范数解 print('X1@b1 == X2@b2:', np.allclose(X1 @ b1, X2 @ b2)) # True print('X1@b1 == X3@b3:', np.allclose(X1 @ b1, X3 @ b3)) # True

b3打印出来会是一组看着挺正常的数，但它的具体取值依赖于库的最小范数约定，换成另一个 SVD 实现可能就不一样了。这就是“能算出来”和“能解释”之间的鸿沟。

条件数这一块也值得单独说。单因子模型里，条件数主要由样本不均衡程度决定。极端一点，假设两组样本量分别是 1 和 9999，设计矩阵 $X=[\mathbf{1}, D_1]$ 的 Gram 矩阵是 $\begin{bmatrix}10000 & 1\ 1 & 1\end{bmatrix}$，特征值约为 10000 和 1，条件数约 $10^4$，而 $X$ 本身的条件数是其平方根约 100。条件数每放大一个数量级，系数的有效数字就掉一位。所以“某个小类别组样本太少导致系数不稳定”，不是玄学，是条件数在说话。

4.2 岭回归为什么能“救”回来，以及救回来的是什么

秩亏的时候加个 $L_2$ 惩罚就能解了，这是常识。但“救回来”的到底是什么，值得掰开说。

岭回归的目标是 $\min_\theta |y-X\theta|^2 + \lambda|\theta|^2$（实际实现里通常不惩罚截距，这里为了公式简洁先写全惩罚，下面代码会处理）。闭式解是：

$$\hat\theta_\lambda = (X^\top X + \lambda I)^{-1}X^\top y$$

即使 $X^\top X$ 奇异，只要 $\lambda>0$，$X^\top X + \lambda I$ 的所有特征值都至少是 $\lambda$，必然正定可逆，解唯一存在。所以岭回归确实把“解不唯一”变成了“解唯一”。

但关键在于：这个唯一解依赖 $\lambda$，也依赖你选的编码方式。$\lambda$ 趋于 0 时，$\hat\theta_\lambda$ 并不收敛到某个“真正的”参数值，而是收敛到最小范数解（准确说是 $X^\top y$ 在 $X^\top X$ 行空间上的投影对应的那个解）。换一套编码，最小范数解就变了。所以岭回归的系数在秩亏场景下不能当作参数估计来解释，它只是“让模型能跑起来”的工程手段。

不过有一件事是稳的：预测值不随编码改变。因为 $X\hat\theta_\lambda$ 落在 $X$ 的列空间里，而列空间对任何满秩重参数化都不变。这一点可以用代码验证：

def ridge_fit(X, y, lam, penalize_intercept=False): p = X.shape[1] P = np.eye(p) * lam if not penalize_intercept: P[0, 0] = 0.0 return np.linalg.solve(X.T @ X + P, X.T @ y) # 换参照组：用 B 当参照 X1b = np.column_stack([np.ones(n), D[['g_A', 'g_C']].values]) for lam in [1e-6, 1e-2, 1.0]: ba = ridge_fit(X1, y, lam) bb = ridge_fit(X1b, y, lam) pa, pb = X1 @ ba, X1b @ bb print(f'lam={lam:g} 参照C系数={np.round(ba,4)} 参照B系数={np.round(bb,4)} ' f'预测最大差={np.abs(pa-pb).max():.2e}')

跑出来的结果会是：随着 $\lambda$ 从 $10^{-6}$ 增大到 1，两套编码的系数差异逐渐缩小（因为惩罚项在“压平”零空间方向上的分量），而预测值的最大差始终在