R语言qr()函数秩判定异常原因及可靠数值秩计算方案问询

基础R的qr()函数默认采用带列主元的QR分解，其内置秩判定逻辑存在固有缺陷：

• 秩判定的阈值计算规则为tol * max(矩阵行数, 矩阵列数) * abs(R[1,1])，即仅以R矩阵左上角第一个对角元（首轮选中的主元）的绝对值作为缩放基准，而非以R矩阵所有对角元的最大绝对值作为基准。
• 本次问题中的badL场景下，矩阵仅1行是真实非零值，其余行均为浮点计算产生的、量级在1e-16~1e-18的机器精度噪声。列主元选择环节受浮点噪声干扰，首轮选中了噪声列作为主元，导致R[1,1]是极小的噪声值，计算出的判定阈值远小于其余噪声对角元的量级，最终将3个噪声主元全部误判为有效主元，输出秩为4的错误结果。
• 调整tol参数无法修正该问题的核心原因是阈值的缩放基准完全错误：无论tol取何值，只要第一主元是噪声值，阈值永远无法匹配真实有效主元的量级。而solve()函数在判定奇异性时是基于全矩阵的条件数计算（通过估计最大/最小奇异值的比值得到），不依赖单主元的缩放，因此能正确识别奇异性。zapsmall()预处理生效的原因是直接将机器精度级的噪声置为0，让列主元选择环节能正确选中真实的非零列作为第一主元，回到正确的判定逻辑上。

• 修正秩判定逻辑，不要直接采信qr()返回的rank字段：完成QR分解后，提取R矩阵的对角元，计算所有对角元的绝对值最大值，以max(dim(x)) * .Machine$double.eps * 对角元绝对值最大值作为判定阈值，统计对角元绝对值大于该阈值的数量作为矩阵的计算秩。该阈值逻辑与条件数判定规则完全一致，目前采用的R矩阵对角元阈值检查思路方向正确，按上述规则做机器精度自适应缩放即可，不需要设置固定阈值。
• 对精度要求高的场景，替换基础QR实现：可以直接调用Matrix包中面向稠密矩阵的QR分解方法，其内置秩判定默认采用R对角元最大值作为缩放基准，不会出现单主元误选导致的秩判定错误；如果对稳定性要求极高、矩阵规模不大，也可以直接用SVD分解做秩判定：将奇异值小于max(dim(x)) * .Machine$double.eps * 最大奇异值的部分截断，结果鲁棒性最强，完全不受主元选择误差影响。
• 优化预处理逻辑，避免粗暴置零：不建议直接用默认参数的zapsmall()做全局预处理，容易误杀量级较小的真实效应。可以在分解前先对矩阵做列归一化（将每列的二范数缩放到1），分解完成后再做逆缩放还原，能大幅降低列量级差异、浮点噪声对主元选择的干扰；如果确实需要清理机器精度噪声，可以将zapsmall()的digits参数设为ceiling(-log10(.Machine$double.eps^0.8))做温和置零，容错性更好。
• 求解线性方程组、计算广义逆时做秩截断：得到正确的秩r后，不要直接用全维度的R矩阵计算，仅提取Q矩阵的前r列、R矩阵的前r行r列的满秩块做计算，得到的广义逆、残差结果不会出现量级极大的异常值，和标准广义逆计算结果一致且运算效率更高，完全适配emmeans包中拆分可估函数、计算自由度和F统计量的需求。

内容的提问来源于stack exchange，提问作者Russ Lenth