dqrsl 拆解：拿着 QR 结果能算出哪 5 种东西**`dqrdc` 把 $Q$ 藏进反射向量里，`dqrsl` 拿

dqrsl 拆解：拿着 QR 结果能算出哪 5 种东西

写在前面：一个函数养活整个最小二乘

上一篇结尾留了个钩子：dqrdc 把 $Q$ 拆成一串 Householder 反射 $P_1 P_2 \cdots P_k$ ，根本不存 $Q$ 这个矩阵，只把每个反射的向量存进 $X$ 的严格下三角、首元编码进 qraux 数组。

于是问题来了：连 $Q$ 都没存，怎么算 $Qy$ ？怎么算 $Q^T y$ ？怎么解最小二乘？

答案就是 dqrsl。它是 dqrdc 的"用户侧"——你给它一份 QR 因子（就是 dqrdc 输出的那份压缩 $X$ 加上 qraux），再给它一个向量 $y$ ，它一次性能给你算出 5 种东西：

qy  = Q·y          左乘正交矩阵
qty = Qᵀ·y         左乘正交矩阵的转置
b   = β̂            最小二乘回归系数
rsd = y − Xβ̂        残差
xb  = Xβ̂ = ŷ        拟合值

一个函数，五种输出。R 的 lm.fit()、SPSS 的 REGRESSION、Python scipy.linalg.lstsq 的底层，全都在调它（或它的等价物）。最小二乘这一整套活，dqrsl 一个函数全包了。今天我们就用真实源码，看它是怎么做到的。

数学原理：反射正向走还是反向走，决定一切

回忆 Householder 反射的定义。对一个单位向量 $u$ ，反射矩阵是

$P = I - 2 u u^T, \quad P^T = P, \quad P^2 = I$

它既对称又正交。dqrdc 把 $X$ 分解成 $X = QR$ ，其中

$Q = P_1 P_2 \cdots P_k, \quad Q^T = P_k \cdots P_2 P_1$

注意这里有两个关键事实，它们决定了 dqrsl 的全部设计：

事实一：因为每个 $P_i$ 都对称， $Q^T = (P_1 P_2 \cdots P_k)^T = P_k \cdots P_2 P_1$ 。 也就是说 $Q$ 是把反射按 $1, 2, \ldots, k$ 的顺序乘起来； $Q^T$ 是把它们反过来乘。

事实二：因为 $P_i^2 = I$ ，把 $P_i$ 作用到一个向量 $z$ 上等于做一次"反射更新"： $P_i z = z - 2 (u_i^T z)\, u_i$ 这就是一行 ddot（算 $u_i^T z$ ）加一行 daxpy（ $z \leftarrow z - 2(u_i^T z) u_i$ ）。 $Q$ 矩阵本身从头到尾没有被构造出来，永远是"作用到向量上"。

于是：

$Q^T y = P_k P_{k-1} \cdots P_1 \, y \quad (\text{正向 } P_1 \to P_k)$ $Q y = P_1 P_2 \cdots P_k \, y \quad (\text{反向 } P_k \to P_1)$

注意一个微妙之处：矩阵乘法右结合， $Q^T y = (P_k \cdots P_1) y$ 在数值上是先作用 $P_1$ 、再 $P_2$ 、……、最后 $P_k$ ；而 $Qy = (P_1 \cdots P_k) y$ 是先作用 $P_k$ 、……、最后 $P_1$ 。作用顺序正好相反。这一点马上会在代码里看到——cqy 分支和 cqty 分支的 for 循环方向一正一反，差之毫厘，结果一个是 $Qy$ 、一个是 $Q^Ty$ 。

至于最小二乘： $X\beta = QR\beta = y$ ，左乘 $Q^T$ 得 $R\beta = Q^T y$ 。而 $R$ 是上三角，回代即可， $O(p^2)$ 。残差和拟合值从 $Q^T y$ 的后 $n-p$ 个和前 $p$ 个分量直接读出来再变换回去。整个最小二乘的计算量被 Householder 结构压到极致。

数学说完了，下面把 qr.c 里的 dqrsl 真实源码一段段贴出来拆。

第一段：job 用十进制位编码选输出

int dqrsl(double *x, int ldx, int n, int k, const double *qraux, const double *y,
          double *qy, double *qty, double *b, double *rsd, double *xb, int job) {
#define X(i,j) x[((i)-1) + ((j)-1)*ldx]
#define QR(k)  qraux[(k)-1]
#define Y(k)   y[(k)-1]
    int info = 0;
    int cqy  = job/10000 != 0;
    int cqty = job%10000 != 0;
    int cb   = (job%1000)/100 != 0;
    int cr   = (job%100)/10 != 0;
    int cxb  = job%10 != 0;
    int ju   = (k < n - 1) ? k : n - 1;            /* min0(k, n-1) */

这是 dqrsl 最有工程美感的地方：一个函数当五个用。

注意 job 是一个普通的 int，但源码注释里写得明明白白：

/10000  QY   ；%10000  QTY   ；(%1000)/100  B   ；(%100)/10  RSD   ；%10  XB

它用十进制的五位来当开关：

位权	含义	解码表达式	想要这个输出，job 至少要包含
10000	QY（ $Qy$ ）	`job/10000 != 0`	10000
1000	QTY（ $Q^Ty$ ）	`job%10000 != 0`	1000
100	B（回归系数 $\hat\beta$ ）	`(job%1000)/100 != 0`	100
10	RSD（残差）	`(job%100)/10 != 0`	10
1	XB（拟合值）	`job%10 != 0`	1

为什么用十进制位而不是二进制位？因为这是 1978 年 Fortran 时代的代码。当年的程序员觉得用十进制位（10000、1000、100……）比二进制位（16、8、4、2、1）更直观、更好从打孔卡片上读。你只要把 job 当成一个五位数 QYQTYBRSDXB 的开关组合，就能一眼看出想要什么——job=11111 全要，job=10100 只要 QY 和 B（实际 LINPACK 文档约定 job 用十进制数字书写，例如 10000 表示只要 QY）。

于是调用者用同一个函数、传不同的 job，按需取输出，不用为每种输出单独写一个函数。这就是为什么我说"一个函数养活整个最小二乘"——它是 1978 年那种内存按 KB 计、函数调用开销敏感的年代，对代码复用的极致压榨。

注意那个 ju = (k < n - 1) ? k : n - 1。这里 k 是 QR 分解的有效秩（实际参与分解的列数），ju 是"真正要用到的反射个数"。因为第 $n$ 个反射 $P_n$ 作用在长度为 1 的子向量上、实际上是恒等（或退化），所以最多用到 $P_{n-1}$ 。这个 min(k, n-1) 在后面的反射循环里到处都用。

第二段：算 $Qy$ ——反向应用反射

if (cqy)  dcopy(n, y, 1, qy, 1);
if (cqty) dcopy(n, y, 1, qty, 1);
/* QY = Q·y ：反向应用各 Householder 反射 */
if (cqy) {
    for (jj = 1; jj <= ju; jj++) {
        j = ju - jj + 1;
        if (QR(j) == 0.0) continue;
        temp = X(j,j); X(j,j) = QR(j);
        t = -ddot(n-j+1, &X(j,j), 1, &qy[j-1], 1) / X(j,j);
        daxpy(n-j+1, t, &X(j,j), 1, &qy[j-1], 1);
        X(j,j) = temp;
    }
}

这段代码有四个细节值得逐字拆。

细节一：先 dcopy。 进入分支前，把 $y$ 复制到 qy（和 qty）缓冲区，后面所有反射都作用在副本上，原始 $y$ 永远不被修改。这是一个干净接口的基本素养。

细节二：循环方向是 j = ju - jj + 1，即从 ju 倒着走到 1。 这就是前面数学推导说的—— $Qy = P_1 P_2 \cdots P_k \, y$ ，矩阵乘法右结合，先作用 $P_k$ ，再作用 $P_{k-1}$ ，……，最后 $P_1$ 。倒着循环就是为了按这个顺序应用。把循环改成正向走，算出来的就不是 $Qy$ 而是 $Q^Ty$ 了。

细节三：临时换对角元再恢复。 这是整段代码最巧妙的工程技巧，注释里专门强调：

temp = X(j,j); X(j,j) = QR(j);    // 应用前：把 qraux 存的首元放回矩阵
t = -ddot(...) / X(j,j);           // 反射更新
daxpy(...);
X(j,j) = temp;                     // 应用后：恢复为 R 的对角元

为什么这么做？回忆 dqrdc 的存储约定：反射向量存在 $X$ 的严格下三角，但它的首元被"挪走"了——对角位置 $X(j,j)$ 存的是 $R$ 的对角元（带符号的 $-NRMXL$ ），而反射向量的真正首元被单独放在 qraux[j] 里。这是为了让一份 $X$ 同时承载 $R$ （显式，在对角和上三角）和反射向量（隐式，在严格下三角，首元外置）。

所以应用反射时，必须先把反射向量完整地拼回来——把 qraux[j] 的首元临时放回 $X(j,j)$ ，让下三角的那段 + 对角位置组成一个完整的 Householder 向量 $u_j$ 。ddot 算 $u_j^T z$ 、daxpy 做 $z \leftarrow z - 2(u_j^T z) u_j / u_j(j)$ （注意除以 $X(j,j)$ ，因为 dqrdc 存的反射向量首元不是 1，需要归一化）。应用完立刻恢复，让 $X(j,j)$ 继续保留 $R$ 的对角元，不破坏后续回代。

这一来一回两行 temp = X(j,j); X(j,j) = QR(j); ... X(j,j) = temp;，是算法工程师对内存的极致节俭——不另开数组存反射向量，复用对角位置，靠临时换入换出实现"一份内存两种用途"。这种 hack 在 1978 年的 Fortran 代码里随处可见，但放到今天看依然精妙。

细节四：反射更新的公式。 把上面三行展开： $u_j^T z \leftarrow \text{ddot}(u_j, z), \quad t = -\frac{u_j^T z}{u_j(j)}, \quad z \leftarrow z + t \cdot u_j$

为什么除以 $u_j(j)$ 又为什么前面有负号？这是 LINPACK 的一种特定归一化约定：dqrdc 存的反射向量 $u_j$ 并没有把首元归一成 1，而是带着一个比例因子；公式 $z \leftarrow z - 2\frac{u_j^T z}{u_j^T u_j} u_j$ 在 LINPACK 的存储约定下被改写成 $z \leftarrow z - \frac{u_j^T z}{u_j(j)} u_j$ （因为 $u_j^T u_j = 2 u_j(j)$ ，正好约掉因子 2）。这是教科书永远不会告诉你的"反射更新公式的生产版本"。

第三段：算 $Q^Ty$ ——正向应用反射

/* QTY = Qᵀ·y ：正向应用（注意方向）*/
if (cqty) {
    for (j = 1; j <= ju; j++) {
        if (QR(j) == 0.0) continue;
        temp = X(j,j); X(j,j) = QR(j);
        t = -ddot(n-j+1, &X(j,j), 1, &qty[j-1], 1) / X(j,j);
        daxpy(n-j+1, t, &X(j,j), 1, &qty[j-1], 1);
        X(j,j) = temp;
    }
}

把这段和上一段对照看——除了循环方向从 j = ju - jj + 1（倒序）变成 j = 1, 2, ..., ju（正序），反射更新的核心三行代码逐字相同：

temp = X(j,j); X(j,j) = QR(j);
t = -ddot(n-j+1, &X(j,j), 1, &qty[j-1], 1) / X(j,j);
daxpy(n-j+1, t, &X(j,j), 1, &qty[j-1], 1);
X(j,j) = temp;

这就是数学原理里说的"反向应用算 $Qy$ 、正向应用算 $Q^Ty$ "。同一组反射向量、同一个更新公式，仅仅换个循环方向，输出就从 $Qy$ 变成 $Q^Ty$ 。这是 Householder 反射相对于 Givens 旋转、相对于显式存 $Q$ 的一个巨大优势——它"作用"和"转置作用"是同一份代码，只差顺序。

理解了这一点你就懂了为什么 dqrsl 能同时算 $Qy$ 和 $Q^Ty$ ：它根本不需要构造 $Q$ 和 $Q^T$ 两个矩阵，它只需要会"按顺序应用反射"。顺序换一下，转置就有了。

这也是为什么 dqrdc 当年敢那么嚣张地"不存 $Q$ "——只要反射向量在， $Q$ 和 $Q^T$ 都是 dqrsl 一次循环的事。

第四段：从 $Q^Ty$ 切出回归系数和残差的"半成品"

if (cb)  dcopy(k, qty, 1, b, 1);
kp1 = k + 1;
if (cxb) dcopy(k, qty, 1, xb, 1);
if (cr && k < n) dcopy(n-k, &qty[kp1-1], 1, &rsd[kp1-1], 1);
if (cxb && kp1 <= n)
    for (i = kp1; i <= n; i++) xb[i-1] = 0.0;
if (cr)
    for (i = 1; i <= k; i++) rsd[i-1] = 0.0;

这段很短，但信息量极大。它利用了一个最小二乘里最优雅的恒等式：

$Q^T y = \begin{pmatrix} R \hat\beta \\ r \end{pmatrix}$

也就是说，把 $y$ 左乘 $Q^T$ 之后得到一个 $n$ 维向量，它的前 $k$ 个分量等于 $R\hat\beta$ （ $R$ 是 $k \times k$ 上三角），后 $n - k$ 个分量等于残差在正交补空间里的表示 $r$ 。

为什么？因为 $X = QR$ ，最小二乘的拟合值 $X\hat\beta = QR\hat\beta$ ，残差 $y - X\hat\beta$ 与 $X$ 的列空间正交（也就是与 $Q$ 的前 $k$ 列正交）。左乘 $Q^T$ 后，前 $k$ 维是 $R\hat\beta$ （拟合值在 $Q$ 的前 $k$ 列上的投影），后 $n-k$ 维是残差在 $Q$ 后 $n-k$ 列上的投影——这正是"残差的本质"。

于是源码做的就是从这个 $Q^Ty$ 里"切两半"：

要回归系数 $\hat\beta$ ：切前 $k$ 个 dcopy(k, qty, 1, b, 1)，下一节回代得到 $\hat\beta$ 。
要残差：切后 $n-k$ 个 dcopy(n-k, &qty[kp1-1], 1, &rsd[kp1-1], 1)，前 $k$ 位补 0，然后下一节反向应用反射，把它从" $Q^T$ 表示"变回原来的坐标系。
要拟合值：切前 $k$ 个到 xb，后 $n-k$ 位补 0，然后反向应用反射得到 $\hat y = X\hat\beta$ 。

注意两个补零的循环：

if (cxb && kp1 <= n)
    for (i = kp1; i <= n; i++) xb[i-1] = 0.0;   // 拟合值后 n-k 位补 0
if (cr)
    for (i = 1; i <= k; i++) rsd[i-1] = 0.0;    // 残差前 k 位补 0

补零不是凑数，是有数学含义的。残差向量在 $Q^T$ 表示下"前 $k$ 维是 0"（因为它正交于列空间），所以补零后再反向作用反射（即左乘 $Q$ ），就恢复出原坐标系下的残差 $y - X\hat\beta$ 。拟合值同理：它在 $Q^T$ 表示下"后 $n-k$ 维是 0"，补零后反向作用反射，就恢复出 $\hat y$ 。

一个 $Q^Ty$ ，切两半、各补零、各反向反射，残差和拟合值都出来了。 这就是为什么残差和拟合值在数值上严格满足 $\hat y + r = y$ ——它们是从同一个 $Q^Ty$ 出发、走对称的两条路算出来的，数值误差天然对消。

第五段：回代求最小二乘解—— $O(p^2)$ 而不是 $O(n^3)$

/* 回代求最小二乘解 b */
if (cb) {
    for (jj = 1; jj <= k; jj++) {
        j = k - jj + 1;
        if (X(j,j) == 0.0) { info = j; break; }
        b[j-1] = b[j-1] / X(j,j);
        if (j != 1) { t = -b[j-1]; daxpy(j-1, t, &X(1,j), 1, b, 1); }
    }
}

这段是教科书意义上的"上三角回代"。此时 b 里存的是 $Q^Ty$ 的前 $k$ 个分量，也就是 $R\hat\beta$ 。要解 $R\hat\beta = c$ 得到 $\hat\beta$ （这里 $c$ 是 $Q^Ty$ 前 $k$ 维），标准做法是从最后一行回代：

$\hat\beta_k = c_k / R(k,k)$ $\hat\beta_j = \frac{c_j - \sum_{i>j} R(j,i) \hat\beta_i}{R(j,j)}, \quad j = k-1, \ldots, 1$

源码就是它的逐行翻译：

外层 j = k - jj + 1 倒序：从第 $k$ 行倒着回代到第 1 行。
b[j-1] = b[j-1] / X(j,j)：除以对角元 $R(j,j)$ 。
daxpy(j-1, t, &X(1,j), 1, b, 1)：把第 $j$ 列前 $j-1$ 行（即 $R(1,j), \ldots, R(j-1,j)$ ）乘上 $-\hat\beta_j$ ，加到 b 的前 $j-1$ 位——这就是 $\sum_{i>j} R(j,i) \hat\beta_i$ 这一项的反向减法。

注意循环范围只有 $k$ （一般 $k = p$ ，列数）。整个回代的计算量是 $O(k^2)$ ，而不是 $O(n^3)$ 。这是 QR 分解相对直接解正规方程的一个隐藏优势：QR 把 $n \times p$ 的最小二乘问题压成一个 $p \times p$ 的上三角系统，剩下的就是 $O(p^2)$ 的回代。当 $n \gg p$ （典型场景：几千个样本、几个变量），这一步几乎免费。

另一个工程细节是 if (X(j,j) == 0.0) { info = j; break; }。如果某个对角元是 0，意味着 $R$ 是奇异的，矩阵 $X$ 不满秩——这是回归里常见的"共线性"问题。dqrsl 不会崩溃，而是把出问题的列号写进 info 返回给调用者。R 的 lm() 在遇到奇异拟合时打印那些被 NA 掉的系数，底层靠的就是这个 info。

第六段：反向反射，恢复残差和拟合值

/* 由残差/拟合反推，应用 Q 得到完整结果 */
if (cr || cxb) {
    for (jj = 1; jj <= ju; jj++) {
        j = ju - jj + 1;
        if (QR(j) == 0.0) continue;
        temp = X(j,j); X(j,j) = QR(j);
        if (cr) { t = -ddot(n-j+1, &X(j,j), 1, &rsd[j-1], 1) / X(j,j);
                  daxpy(n-j+1, t, &X(j,j), 1, &rsd[j-1], 1); }
        if (cxb){ t = -ddot(n-j+1, &X(j,j), 1, &xb[j-1], 1) / X(j,j);
                  daxpy(n-j+1, t, &X(j,j), 1, &xb[j-1], 1); }
        X(j,j) = temp;
    }
}

这一段和前面的 $Qy$ 分支几乎一模一样——倒序循环、临时换对角元、反射更新、恢复。区别只是作用对象从 qy 变成了 rsd 和 xb（两个并行做）。

数学含义是：残差和拟合值在前面被切成了 $Q^T$ 表示下的"半成品"（一半补零），现在反向应用反射、也就是左乘 $Q$ ，把它们变回原来的坐标系。

注意循环倒序——j = ju - jj + 1。回忆前面的规律：倒序应用反射算的是 $Q \cdot (\text{向量})$ 。残差半成品是 $Q^T(y - X\hat\beta) = (0, \ldots, 0, r_{k+1}, \ldots, r_n)^T$ ，再左乘 $Q$ 就得到原坐标系下的残差 $y - X\hat\beta$ 。拟合值半成品同理，左乘 $Q$ 得到 $X\hat\beta$ 。

把这一段和第二段（算 $Qy$ ）对照看：两段代码框架完全相同，只是作用对象不同。这是 dqrsl 代码复用的另一个体现——同一个"反向应用反射"的工具，既用来从 $y$ 算 $Qy$ ，也用来从半成品算出残差和拟合值。

至此，五种输出全部算完。回归系数走回代，残差/拟合值走"切半 + 反向反射"， $Qy$ 和 $Q^Ty$ 走"反射的正反两种顺序"。一份 QR 因子，全部拿下。

第七段：退化情形与返回值

} else {                                        /* ju == 0 的退化情形 */
    if (cqy)  qy[0]  = Y(1);
    if (cqty) qty[0] = Y(1);
    if (cxb)  xb[0]  = Y(1);
    if (cb) {
        if (X(1,1) == 0.0) info = 1; else b[0] = Y(1) / X(1,1);
    }
    if (cr) rsd[0] = 0.0;
}
#undef X
#undef QR
#undef Y
    return info;
}

ju == 0 是退化情形：没有反射需要应用（ $k = 0$ 或者 $n = 1$ 的边界）。这时 $Q = I$ ，所以 $Qy = Q^Ty = y$ ，残差为 0、拟合值就是 $y$ 。回归系数直接由 $R(1,1)\hat\beta_1 = y_1$ 解出（如果 $R(1,1) = 0$ ，记 info = 1 报告秩亏）。

最后的 return info 是整个函数的唯一对外信号——0 表示成功，非 0 表示第 info 列对角元为零（即 $X$ 秩亏）。这就是调用者判断"这次最小二乘是否可信"的唯一依据。

教科书不讲、但工程上至关重要的 5 个点

把整段源码读完，我们能总结出 5 个教科书几乎不提、但真实代码里全是心血的工程细节。

1. 同一份 QR 结果算 5 种东西，用十进制位编码做开关。 一个 dqrsl 函数，传不同的 job（10000 / 1000 / 100 / 10 / 1），按需输出 QY / QTY / B / RSD / XB 中的任意组合。这种"一个函数当五个用"的设计是 1978 年 Fortran 时代对代码复用的极致压榨，今天的 R lm.fit 和 SciPy 的 lstsq 底层依然这么调。

2. 反射的正向/反向应用决定算 $Qy$ 还是 $Q^Ty$ 。 同一组 Householder 向量、同一个反射更新公式，仅仅换个循环方向（j = 1..ju 还是 j = ju..1），结果就从 $Qy$ 变成 $Q^Ty$ 。这就是 Householder 反射不需要单独构造 $Q$ 矩阵的根本原因——它的"作用"和"转置作用"是同一份代码。

3. "临时换对角元再恢复"的技巧。 反射向量首元被外置到 qraux，应用反射前必须把它临时放回 $X(j,j)$ 、应用完立刻恢复，让 $X$ 既承载 $R$ （显式）又承载反射向量（隐式）。这种"一份内存两种用途"的 hack 是内存按 KB 计年代留下的智慧，今天看依然精妙。

4. 回归系数只需 $R$ 的前 $k$ 行回代， $O(p^2)$ 而非 $O(n^3)$ 。 QR 分解把 $n \times p$ 的最小二乘压成 $p \times p$ 上三角系统，剩下的就是 $O(p^2)$ 回代。当 $n \gg p$ （典型统计场景），这一步几乎免费。

5. 残差和拟合值从同一个 $Q^Ty$ 切出来，数值上严格互补。 一个 $Q^Ty$ 切成前 $k$ 维（拟合值半成品）和后 $n-k$ 维（残差半成品），各自补零后反向反射。因为是从同一个向量出发、走对称的两条路，残差和拟合值天然满足 $\hat y + r = y$ ，数值误差对消。这就是为什么 SPSS/R 报告的残差和拟合值加起来分毫不差地等于 $y$ 。

把整条链串起来：dqrsl 支撑了什么

读懂 dqrsl 这 80 行，你就读懂了统计软件回归功能的内核：

dqrsl（应用 QR 因子）
   ├── 回归系数 β̂：R 回代（O(p²)）
   ├── 拟合值 ŷ = Xβ̂：QTy 切前 k 维 + 反射
   ├── 残差 y − Xβ̂：QTy 切后 n−k 维 + 反射
   ├── Qy / Qᵀy：反射正反两种顺序应用
   └── 由此派生的统计量：
        ├── 杠杆值（hat matrix 对角元，由 Q 的前 k 列范数平方）
        ├── 标准化残差、学生化残差
        ├── ANOVA 平方和（残差向量的范数²）
        ├── 置信区间 / 预测区间
        └── Cook 距离、DFBETAS 等影响度量

R 的 lm.fit() 把 job 设成 11111（全要），一次调用拿到 $\hat\beta$ 、 $\hat y$ 、残差、QR 分解本身；SPSS 的 REGRESSION 命令、UNIANOVA 命令底层都是这条路径。最小二乘这一整套活，dqrsl 一个函数全包了。

更妙的是，这个函数和上一篇的 dqrdc 构成了一个完美的二段式管线：dqrdc 负责"分解、压缩、不存 $Q$ "，dqrsl 负责"按需解压、按需计算"。两个函数加起来不到 200 行 C 代码，构成了过去 40 多年所有主流统计软件回归功能的算法核心。1978 年 LINPACK 的工程师用最小的代码量、最巧的存储约定、最精的反射应用，把最小二乘这件事彻底解决——今天的我们打开任何一本统计计算教材，看到的依然是这套方案的变体。

总结：教科书不讲的 5 个点

一份 QR 算 5 种东西：QY / QTY / B / RSD / XB，通过 job 的十进制位编码选择，一个函数当五个用。
反射的方向决定一切：正向应用算 $Q^Ty$ ，反向应用算 $Qy$ ，同一组向量、同一份公式，只差循环顺序。
临时换对角元再恢复：反射向量首元外置到 qraux，应用时临时放回、应用完恢复，一份内存两种用途。
回归系数 $O(p^2)$ 回代：QR 把 $n \times p$ 压成 $p \times p$ 上三角， $n \gg p$ 时这一步几乎免费。
残差与拟合值互补：从同一个 $Q^Ty$ 切出，数值上严格满足 $\hat y + r = y$ 。

如果只记一句话，那就是：dqrdc 把 $Q$ 藏进反射向量里，dqrsl 拿着这份压缩包按需解压——整个最小二乘只需要 80 行 C 代码就能跑通。 这不是抽象的数学课，而是 1978 年留下来、至今仍在每一台跑统计软件的电脑上日复一日运行的具体代码。

下一篇，我们会继续在 QR 家族里转，看看列主元 QR 是怎么在不增加额外存储的前提下、顺手把矩阵的秩和共线性问题一起解决的。那又是一段教科书糊弄过去、代码里精妙绝伦的故事。

dqrsl 拆解：拿着 QR 结果能算出哪 5 种东西