列主元 QR 与数值秩证书

形式陈述

给定实稠密矩阵 $A\in\mathbb R^{m\times n}$ ， $m\ge n$ ，列主元 QR 在QR 分解中加入列置换，得到

AP=QR,\qquad Q^{\mathsf T}Q=I_m,\qquad R=\begin{pmatrix}R_{\mathrm{top}}\\0\end{pmatrix},

其中 $P$ 是 $n\times n$ 置换矩阵， $R_{\mathrm{top}}$ 是 $n\times n$ 上三角矩阵。本页用完整 $Q$ 说明分块；实际保存反射向量即可，无须形成 $m\times m$ 矩阵。若置换向量为 $p$ ，约定 $P=(e_{p_1},\ldots,e_{p_n})$ ，因此 $AP$ 的第 $j$ 列是原来的 $a_{p_j}$ 。

算法每步选择尚未被已有方向解释的最长列。具体地，令工作矩阵 $B=A$ ， $p=(1,\ldots,n)$ 。第 $j$ 步执行：

在 $\ell=j,\ldots,n$ 中选择使 $\|B_{j:m,\ell}\|_2$ 最大的 $\ell$ ；并列时任取。
交换 $B$ 的第 $j,\ell$ 整列和对应的 $p_j,p_\ell$ 。前 $j-1$ 行已经存好的系数也必须跟随列交换。
对 $x=B_{j:m,j}$ 构造Householder 反射。若 $x\ne0$ ，取 $\alpha=-\operatorname{sign}(x_1)\|x\|_2$ （ $x_1=0$ 时取正号）， $v=x-\alpha e_1$ ， $H=I-2vv^{\mathsf T}/(v^{\mathsf T}v)$ ；将 $H$ 左乘当前尾块 $B_{j:m,j:n}$ 。若 $x=0$ ，本步用恒等变换。
保存反射向量，将主元列对角线下方视为已消去；更新下一步各列的尾范数。精确算术中可从旧范数平方减去新产生的第 $j$ 行元素平方；浮点中若相减严重抵消，应从实际尾列重新计算。

这里参与选择的是当前尾部，不是输入时的原始列范数。反射已将先前选出的方向搬到前几行，删除这些行后留下的长度，才是到已选列空间的正交距离。实现还应以稳定范数例程处理尺度，检查非有限输入；上述反射公式说明数学操作，直接平方求和可能溢出。

算法输出隐式 $Q$ 、 $R$ 和 $p$ 。若还要报告数值秩，则输入必须另含绝对容差 $\tau\ge0$ ，或由相对容差与矩阵尺度换算出的 $\tau$ 。定义

r_\tau(A)=\#\{i:\sigma_i(A)>\tau\}.

对候选 $1\le k<n$ ，把三角因子分块为

R=\begin{pmatrix}R_{11}&R_{12}\\0&R_{22}\end{pmatrix}, \qquad R_{11}\in\mathbb R^{k\times k}.

这里 $R_{22}$ 包括完整 $R$ 的底部零行。保留前 $k$ 个正交方向构造

A_k=Q_{:,1:k}\begin{pmatrix}R_{11}&R_{12}\end{pmatrix}P^{\mathsf T}.

由谱范数的正交不变性， $\|A-A_k\|_2=\|R_{22}\|_2$ ，并且有两个可以共同使用的界：

\sigma_k(A)\ge\sigma_{\min}(R_{11}),\qquad \sigma_{k+1}(A)\le\|R_{22}\|_2.

因此，若

\boxed{\ \|R_{22}\|_2\le\tau<\sigma_{\min}(R_{11})\ },

便得到 $r_\tau(A)=k$ 的证书。只检验右下角很小，只能说明存在一个好的低秩近似，尚不能说明保留下来的 $k$ 个方向都高于阈值。

两个不等式为什么成立

置换后的前 $k$ 列组成 $C=Q_{:,1:k}R_{11}$ ，故 $\sigma_k(C)=\sigma_{\min}(R_{11})$ 。把其余列记作 $D$ ，则 $AA^{\mathsf T}=CC^{\mathsf T}+DD^{\mathsf T}$ 。后一项半正定，特征值的极小极大原理说明第 $k$ 大特征值不会减小，从而给出第一个界。这里使用SVD将奇异值平方识别为 $AA^{\mathsf T}$ 的特征值。

第二个界来自同一 SVD 工具的最佳低秩逼近结论：所有秩至多 $k$ 的矩阵与 $A$ 的谱范数距离，最小值是 $\sigma_{k+1}(A)$ 。上述 $A_k$ 是其中一个候选，故 $\sigma_{k+1}(A)\le\|A-A_k\|_2=\|R_{22}\|_2$ 。两侧同时跨过 $\tau$ ，才排除了第 $k$ 个方向太弱和第 $k+1$ 个方向太强这两种可能。

直觉

每加入一列，已有列空间就多解释一个方向。下一步应选离这个空间最远的剩余列：一根原本很长、却几乎完全落在旧空间里的列，新增信息可能很少。尾范数正是扣除已有解释之后的大小；列置换使这些较新的方向尽早进入 $R_{11}$ 。

分块有两个不同职责。 $R_{22}$ 衡量忽略后续方向会损失多少， $R_{11}$ 衡量保留方向能否可靠地分开。小尾块并不自动带来良好的保留块。按照条件数与扰动的观点，求解系数还应报告 $\kappa_2(R_{11})$ ，而不能仅以“秩是 $k$ ”代替敏感度分析。

例子与边界

从主元到证书的完整例子

取

A=\begin{pmatrix}1&0&2\\1&2&2\\\varepsilon&0&0\end{pmatrix}, \qquad 0<\varepsilon<\sqrt2.

原始三列的长度依次为 $\sqrt{2+\varepsilon^2},2,\sqrt8$ ，故先选 $a_3$ 。可将第一正交方向取为 $q_1=(1,1,0)^{\mathsf T}/\sqrt2$ 。扣除在 $q_1$ 上的投影后， $a_1$ 剩下 $\varepsilon e_3$ ，长度为 $\varepsilon$ ； $a_2$ 剩下 $(-1,1,0)^{\mathsf T}$ ，长度为 $\sqrt2$ 。所以第二步选 $a_2$ ，最后才选 $a_1$ ，得到 $p=(3,2,1)$ 。

调整反射产生的列符号，使三角对角元为正，可以写出

Q=\begin{pmatrix} 1/\sqrt2&-1/\sqrt2&0\\ 1/\sqrt2&1/\sqrt2&0\\ 0&0&1 \end{pmatrix},\qquad R=\begin{pmatrix} 2\sqrt2&\sqrt2&\sqrt2\\ 0&\sqrt2&0\\ 0&0&\varepsilon \end{pmatrix}.

逐列相乘给出 $QR=(a_3,a_2,a_1)=AP$ 。在 $k=2$ 处分块：

R_{11}=\sqrt2\begin{pmatrix}2&1\\0&1\end{pmatrix},\quad R_{12}=\begin{pmatrix}\sqrt2\\0\end{pmatrix},\quad R_{22}=(\varepsilon).

$R_{11}^{\mathsf T}R_{11}=\begin{pmatrix}8&4\\4&4\end{pmatrix}$ 的特征值为 $6\pm2\sqrt5$ ，因此

\sigma_{\min}(R_{11})=\sqrt{6-2\sqrt5}\approx1.23607,\qquad \kappa_2(R_{11})=\frac{3+\sqrt5}{2}\approx2.61803.

取 $\varepsilon=10^{-6}$ ， $\tau=10^{-4}$ ，就有

\sigma_3(A)\le10^{-6}<10^{-4}<1.23607\le\sigma_2(A).

所以 $r_\tau(A)=2$ ，而 $\det A=-4\varepsilon\ne0$ ，精确秩仍为 $3$ 。若希望按输入尺度报告，同一阈值可写成 $\tau/\|A\|_F=10^{-4}/\sqrt{14+\varepsilon^2}$ ；换一个容差是在换有效秩问题。

还可以明确指出哪一列近似多余。对 $C=(a_3,a_2)$ ，用三角求解得到 $R_{11}T=R_{12}$ 的解 $T=(1/2,0)^{\mathsf T}$ ，于是

a_1=CT+\varepsilon e_3=\tfrac12a_3+\varepsilon e_3.

这同时给出所选列、其余列的表示系数和大小恰为 $\varepsilon$ 的剩余误差。

主元规则能说明多少

精确算术下，贪心选择使 $|r_{11}|\ge\cdots\ge|r_{nn}|$ ，但这些对角元不是奇异值。普通 CPQR 在特制矩阵上可能不能充分揭示谱间隙；Gu–Eisenstat 的 Kahan 型反例说明，不能从该贪心规则无条件推出强的秩揭示保证。即使每个剩余列的尾范数至多为 $\delta$ ，直接得到的也只是

\|R_{22}\|_2\le\|R_{22}\|_F\le\sqrt{n-k}\,\delta,

并非 $\|R_{22}\|_2\le\delta$ 。当证书两侧未能夹住阈值时，结论是本次分块不足以确认该数值秩，可进一步计算相关奇异值或改用更强的列选择。

浮点输出还需检查 $\|AP-\widehat Q\widehat R\|$ 和 $\|\widehat Q^{\mathsf T}\widehat Q-I\|$ 。上述精确不等式不能直接套在不完全正交的 $\widehat Q$ 上作为严格证书。若已建立 $A+E=\widetilde Q\widehat R P^{\mathsf T}$ 、 $\widetilde Q$ 正交及 $\|E\|_2\le\eta$ 的有效后向误差界，则奇异值扰动界给出

\sigma_k(A)\ge\sigma_{\min}(\widehat R_{11})-\eta, \qquad \sigma_{k+1}(A)\le\|\widehat R_{22}\|_2+\eta.

严格浮点证书还要为两个块的范数和最小奇异值计算提供可靠上下界；通常的数值估计和小重构残差是诊断，不能自行补齐这些保证。列缩放也会改变奇异值和 $r_\tau$ ，应明确秩判断针对原矩阵还是缩放后的模型。

推论与应用

成本与截断

维护尾范数的实稠密 CPQR，执行前 $k$ 步的主项约为

4mnk-2(m+n)k^2+\frac43k^3

次浮点运算；令 $k=n$ 得到熟悉的 $2mn^2-2n^3/3$ 。该式不包括显式生成 $Q$ 、异常情况下的范数重算和额外诊断成本。可以在部分分解后检查尾块，但此时它只是尚未三角化的矩形块；误差分块公式依然成立。完整分解例程如 LAPACK GEQP3 返回因子与置换，选 $k$ 和解释阈值仍是调用方的工作。

从列选择到求解的边界

若 $R_{11}$ 可逆，解 $R_{11}T=R_{12}$ 后，置换后的矩阵有近似列表示

AP\approx C\begin{pmatrix}I_k&T\end{pmatrix},\qquad \left\|AP-C\begin{pmatrix}I_k&T\end{pmatrix}\right\|_2=\|R_{22}\|_2.

计算时应求解三角系统而非形成逆矩阵。较大的 $T$ 会使列表示对系数误差更敏感。强秩揭示 QR 会进行额外换列；Gu–Eisenstat 的结果在参数 $f>1$ 下同时控制 $|T_{ij}|\le f$ ，并以 $g=\sqrt{1+f^2k(n-k)}$ 给出 $\sigma_i(R_{11})\ge\sigma_i(A)/g$ 与 $\sigma_j(R_{22})\le g\sigma_{k+j}(A)$ 。这是比普通 CPQR 更强的算法保证，但不能消除输入本身很小的奇异值。

对最小二乘，满列秩时可解 $Ry=Q^{\mathsf T}b$ 的前 $n$ 个方程，再由 $x=Py$ 撤销列置换。秩为 $k<n$ 时，令其余变量为零并解 $R_{11}y_1=c_1$ ，只得到一个基本解：它对精确秩亏系统是最小二乘解，对截断模型是该模型的最小二乘解，一般不具有最小范数。

最小例子是 $A=\begin{pmatrix}1&1\\0&0\end{pmatrix}$ 、 $b=(1,0)^{\mathsf T}$ 。基本解 $(1,0)^{\mathsf T}$ 与 $(1/2,1/2)^{\mathsf T}$ 都给零残差，后者范数更小。需要最小范数时，应进一步做完全正交分解或使用 SVD；列主元选择本身没有完成这一步。

参考资料

Peter Businger and Gene H. Golub, Linear Least Squares Solutions by Householder Transformations, Numerische Mathematik 7, 1965, pp. 269–276，尤其 §§1、4：列主元选择、范数更新与置换记录。
Ming Gu and Stanley C. Eisenstat, Efficient Algorithms for Computing a Strong Rank-Revealing QR Factorization, SIAM Journal on Scientific Computing 17(4), 1996, pp. 848–869，§2 的普通 CPQR 与反例、§3 Theorem 3.2 的强秩揭示界。
LAPACK Users’ Guide, 3rd ed., 1999, QR Factorization with Column Pivoting：分解与秩判断的职责，以及基本解的最小范数边界。

形式陈述 ​

两个不等式为什么成立 ​

从主元到证书的完整例子 ​

主元规则能说明多少 ​

成本与截断 ​

从列选择到求解的边界 ​

反向链接