1. 为什么今天还要手撕 QR 分解一个被低估的“数值稳定器”你有没有遇到过这样的情况用 R 或 Python 做线性回归数据一多、变量一杂模型系数突然飘得离谱——明明特征之间只是轻微相关lm()或LinearRegression()却报出一堆极大正负值甚至Inf或NaN或者在做主成分分析前对设计矩阵做 SVD结果前几个奇异值还算稳后面全变成接近零的抖动噪声又或者在 PySpark 上跑大规模最小二乘任务反复失败日志里只有一行模糊的LAPACK error这些不是玄学是浮点数在向你敲门。而 QR 分解就是那个站在浮点误差风暴最前沿、默默扛下所有冲击的“数值稳定器”。它不炫技不抢镜但几乎所有现代统计计算底层都在用它——R 的lm()默认走 QR 路径NumPy 的lstsq底层调用 LAPACK 的dgeqrfSpark MLlib 的LinearRegression在分布式 QR 上做了大量工程优化。它不是教科书里一个冷冰冰的定理而是你每次fit()背后那根最结实的承重梁。我带过三届数据科学训练营每届都有学员卡在“为什么我的回归系数和别人差十倍”这个问题上。翻来覆去查数据、查缺失值、查标准化最后发现根源是他们用的是(X.T X).I X.T y这种“教科书式”解法。这个公式在数学上完全正确但在计算机里X.T X会把原本微小的条件数放大平方倍。举个具体例子假设原始设计矩阵 X 的条件数是 100已经算中等病态X.T X的条件数就飙升到 10,000。此时哪怕输入数据只有 1e-8 的舍入误差解出来的系数也可能偏差 100%。而 QR 分解绕开了X.T X这个“误差放大器”直接在原始空间里操作把数值稳定性从“听天由命”拉回到“可预期、可控制”的水平。这正是它不可替代的核心价值它不改变问题的数学本质却彻底改变了问题在有限精度机器上的求解路径。它不是为了解出“理论最优解”而是为了让你在真实世界的数据、真实的硬件、真实的内存限制下拿到那个“足够好、足够稳、足够可信”的解。所以当你看到“QR Decomposition in Machine Learning”这个标题时请别把它当成一个待背诵的算法名词。它是一把钥匙一把打开高维数值计算可靠性大门的钥匙它是一面镜子照见我们日常所用的lm()、fit()、solve()等函数背后最朴素也最坚韧的工程哲学——在不确定的世界里用确定的数学结构去对抗不确定的误差。2. 核心原理拆解Q 和 R 到底在“分解”什么矩阵分解的本质是给一个复杂的、难以直接处理的数学对象找到一组更简单、更有规律、更易操控的“积木”。QR 分解的精妙之处在于它选的这两块积木——Q 和 R——恰好拥有截然不同、却又完美互补的数学性质。理解它们各自的“性格”比死记A QR这个等式重要十倍。2.1 Q 矩阵那个“不扭曲世界”的守门人Q 是一个正交矩阵。这个词听起来很抽象但它的物理意义极其直观它代表的是一种纯粹的旋转或反射不拉伸、不压缩、不剪切。想象你有一张印着坐标网格的透明胶片上面画着几个向量。现在你把这个胶片在三维空间里任意旋转比如绕 X 轴转 30 度再绕 Y 轴转 45 度或者对着某一面镜子翻转一下。做完这些动作后胶片上的所有向量长度没变任意两个向量之间的夹角也没变它们的点积内积原封不动地保留了下来。这就是正交矩阵 Q 的全部秘密。它的数学定义Q^T Q I翻译成大白话就是“我把一个向量 v 乘以 Q 得到新向量 u再把 u 乘以 Q 的转置就又变回了原来的 v。” 这说明 Q 的逆矩阵就是它自己的转置求逆这件事变得异常廉价——不需要任何复杂的计算只要把矩阵翻过来就行。更重要的是||Qv|| ||v||长度不变和u^T v (Qu)^T (Qv)点积不变这两个性质保证了所有基于向量长度和角度的几何关系在经过 Q 变换后都毫发无损。在机器学习里这个性质意味着什么意味着当我们用 Q 去“预处理”你的设计矩阵 X 时即计算Q^T X我们没有丢失任何关于数据“形状”的信息。X 的列向量之间的相关性、它们的相对距离、它们构成的子空间……所有这些核心统计结构都被 Q 完美地、无损地保存了下来。它就像一个高保真的数字转换器把原始数据从一个可能“嘈杂”的坐标系投射到一个全新的、彼此正交的坐标系里为后续的计算铺平了道路。2.2 R 矩阵那个“秩序井然”的收纳盒如果说 Q 是一个优雅的舞者那么 R 就是一个极度理性的收纳师。它是一个上三角矩阵这意味着它所有的“故事”都写在对角线及其上方。对角线以下的所有元素统统为零。这种结构看似简单实则蕴含着巨大的计算红利。上三角矩阵最大的好处是求解线性方程组Rβ c变得像剥洋葱一样容易。你从最后一行开始R_{nn} * β_n c_n直接就能解出β_n c_n / R_{nn}然后代入倒数第二行R_{n-1,n-1} * β_{n-1} R_{n-1,n} * β_n c_{n-1}因为β_n已知所以β_{n-1}也能立刻算出来如此往复一路向上直到解出β_1。这个过程叫做“回代法”Back Substitution其时间复杂度仅为 O(n²)远低于通用矩阵求逆的 O(n³)。在 QR 分解中R 承载着原始矩阵 A 的全部“尺度”和“方向”信息但它把这些信息以一种高度结构化的方式打包好了。R的对角线元素R_{ii}本质上就是第 i 个正交基向量q_i在原始数据 A 的列空间中所“占据”的“分量大小”。而R的上三角部分则编码了各个正交基向量之间如何“协同”地重构出原始数据。因此R 不是一个随意的中间产物它是整个分解过程的“成果报告”清晰地告诉你为了重建 A你需要多少份q_1多少份q_2以及q_1和q_2之间需要怎样的组合比例。2.3 为什么是 Q 和 R 的组合——一场关于“解耦”的精密手术现在把 Q 和 R 放在一起看A QR。这个等式揭示了一场精妙的“解耦”手术。左边的 A是一个混杂了旋转、缩放、相关性的“混沌体”。右边的 QR则把这个混沌体干净利落地拆成了两个独立的部分Q 部分负责“姿态”它告诉你原始数据所在的那个子空间是以什么样的“朝向”orthonormal basis存在的。R 部分负责“内容”它告诉你在这个已经梳理好的、彼此正交的“朝向”下原始数据的具体“数值构成”是什么。这种解耦带来了革命性的计算优势。当我们要求解最小二乘问题min ||Ax - b||²时传统方法要面对A^T A x A^T b这个病态方程。而用 QR 分解问题瞬间被转化min ||QRx - b||² min ||Rx - Q^T b||²。因为 Q 是正交的||QRx - b||² ||Rx - Q^T b||²这个等式恒成立。于是我们只需要解一个上三角系统Rx Q^T b而Q^T b的计算不过是一系列高效的向量点积。整个过程避开了A^T A这个潜在的“雷区”将一个可能数值不稳定的全局优化问题降维成一个绝对稳定的局部求解问题。这就像修理一台精密钟表。传统方法是把整个机芯拆开试图在混乱的齿轮堆里直接调整某个游丝的张力而 QR 方法则是先用一套标准的校准工具Q把所有齿轮的基准位置都重新标定好然后再针对那个已经孤立出来的、结构清晰的游丝R进行精准微调。前者费力且易出错后者高效且可靠。3. 三种主流实现方式深度对比Gram-Schmidt、Householder、GivensQR 分解不是只有一个“标准答案”而是一套“工具箱”。不同的工具适用于不同的场景它们在精度、速度、内存占用和实现难度上各有千秋。选择哪个不是看谁更“高级”而是看你的数据长什么样、你的机器有多强、你的代码要跑多久。3.1 Gram-Schmidt 过程最直观的“手工课”也是最脆弱的“教学模型”Gram-Schmidt格拉姆-施密特过程是理解 QR 分解的“启蒙老师”。它的思想无比朴素就像搭积木一样一个向量一个向量地构建正交基。第一步取 A 的第一列a₁把它单位化得到q₁ a₁ / ||a₁||。这是我们的第一个正交基向量。第二步取 A 的第二列a₂计算它在q₁方向上的投影proj_q₁(a₂) (a₂·q₁) q₁然后用a₂减去这个投影得到一个与q₁垂直的新向量v₂ a₂ - proj_q₁(a₂)再把它单位化得到q₂ v₂ / ||v₂||。第三步对a₃依次减去它在q₁和q₂方向上的投影得到v₃再单位化得q₃。如此循环直到处理完所有列。这个过程在 R 代码里写出来逻辑清晰得像一篇散文。它完美地诠释了“正交化”和“单位化”的每一个步骤是建立直觉的绝佳途径。提示我在第一次手写 Gram-Schmidt 时犯了一个经典错误在计算v₂后我直接用v₂去计算v₃的投影而不是用已经单位化的q₂。这导致q₃并不真正与q₂正交。正确的做法永远是用当前已有的、单位化的q_i去投影下一个向量。然而Gram-Schmidt 的致命弱点在于数值不稳定性。在浮点运算中当a₂和q₁非常接近即高度相关时a₂ - proj_q₁(a₂)这个减法操作会产生严重的“有效数字抵消”。两个几乎相等的大数相减结果的有效位数会急剧减少导致v₂的精度严重受损进而污染后续所有q_i的计算。这就是为什么你在原文中看到自己手写的 Gram-Schmidt 版本得到的 R 矩阵对角线下方并非严格的零而是充满了1e-15、1e-16这样的“幽灵数字”而 R 内置的qr()函数给出的 R 矩阵则干净利落。这不是你代码有 bug而是算法本身的“体质”决定的。3.2 Householder 反射工业级的“稳定之王”LAPACK 的默认选择如果你把 Gram-Schmidt 比作一位耐心的手工匠人那么 Householder 反射就是一台高精度的数控机床。它的核心思想是我不需要一步步“削”出正交基我只需要用一系列完美的“镜面反射”一次性把矩阵 A 的下三角部分“砸”成零。一个 Householder 矩阵H I - 2uu^T其中u是一个单位向量代表一次关于超平面的反射。它的神奇之处在于对于任何一个向量x你总能找到一个合适的u使得Hx的结果是x关于某个特定超平面的镜像并且这个镜像可以被精确地设计成——让x的后n-k个分量全部变为零。Householder 的实现流程是迭代的对矩阵 A 的第一列a₁构造一个 Householder 矩阵H₁使得H₁a₁的结果是一个只有第一个元素非零的向量[α, 0, 0, ..., 0]^T。将H₁作用于整个矩阵 A得到H₁A此时H₁A的第一列下方全为零。接下来对H₁A的右下角子矩阵去掉第一行第一列重复这个过程构造H₂并作用于H₁A得到H₂H₁A此时前两列的下方都为零。如此继续直到所有下三角元素都被“反射”为零。Householder 的优势是压倒性的极高的数值稳定性它避免了 Gram-Schmidt 中危险的减法操作所有计算都基于向量的外积和矩阵乘法对舍入误差天然免疫。计算效率高虽然每次反射需要计算一个完整的矩阵H但实际应用中H并不显式存储而是只存储向量u所有运算都围绕u展开大大节省了内存和计算量。易于并行化每个 Householder 反射的计算都是独立的非常适合现代 CPU 和 GPU 的并行架构。这也是为什么 R 的qr()函数、Python 的numpy.linalg.qr()其底层都调用 LAPACK 的dgeqrf例程——它正是基于 Householder 反射。它不是“理论上好”而是经过数十年、数十亿次生产环境验证的“实践之王”。3.3 Givens 旋转稀疏矩阵的“外科医生”精准打击每一处病灶Givens 旋转是三种方法里最“精细”的一位。它的武器不是一次覆盖整列的“反射”而是针对单个元素的“旋转”。一个 Givens 矩阵G(i, j, θ)是一个单位矩阵只是在第i行第i列、第i行第j列、第j行第i列、第j行第j列这四个位置上填入了cosθ和sinθ形成一个二维平面内的旋转。它的操作目标非常明确给定矩阵 A 中的两个元素a_pq和a_rq同一列不同行我能否找到一个旋转角度θ使得G * A的结果中a_rq这个位置的值变成零答案是肯定的。通过简单的三角函数计算tan(2θ) 2*a_pq*a_rq / (a_pq² - a_rq²)就能解出θ。Givens 的优势在于极致的稀疏性友好。因为它每次只修改矩阵的两行所以当你的原始矩阵 A 本身就很稀疏比如推荐系统中的用户-物品交互矩阵时Givens 旋转能最大限度地保持这种稀疏性。它不会像 Householder 那样一次反射就把一整列都“染色”从而产生大量新的非零元fill-in。在处理超大规模稀疏矩阵时Givens 的内存占用和计算开销往往能比 Householder 低一个数量级。注意Givens 旋转的缺点是为了将一列下方的所有元素清零它需要n-1次独立的旋转操作而 Householder 只需要一次。所以对于稠密矩阵Givens 的总计算量更大。它的舞台永远属于那些“大而稀”的数据。4. 实操详解从零开始在 R 中完成一次完整的 QR 分解纸上得来终觉浅绝知此事要躬行。下面我将带你用 R 语言亲手完成一次从数据准备、分解、验证到应用的全流程实操。我们不再依赖qr()的黑盒而是用最基础的向量和矩阵运算一步步构建出属于你自己的 QR 分解。4.1 数据准备与探索理解你的“病人”首先让我们加载并审视原文中提供的那个 32x3 行的matrix_A。这段数据看起来像是汽车数据集mpg、weight、hp但它的核心价值不在于业务含义而在于它是一个典型的、具有现实挑战性的“小而病态”样本。# 创建原始矩阵 matrix_A - matrix(c( 21.0, 21.0, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8, 16.4, 17.3, 15.2, 10.4, 10.4, 14.7, 32.4, 30.4, 33.9, 21.5, 15.5, 15.2, 13.3, 19.2, 27.3, 26.0, 30.4, 15.8, 19.7, 15.0, 21.4, 2.62, 2.875, 2.32, 3.215, 3.44, 3.46, 3.57, 3.19, 3.15, 3.44, 3.44, 4.07, 3.73, 3.78, 5.25, 5.424, 5.345, 2.2, 1.615, 1.835, 2.465, 3.52, 3.435, 3.84, 3.845, 1.935, 2.14, 1.513, 3.17, 2.77, 3.57, 2.78, 110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, 180, 180, 205, 215, 230, 66, 52, 65, 97, 150, 150, 245, 175, 66, 91, 113, 264, 175, 335, 109 ), nrow 32, ncol 3, byrow FALSE) colnames(matrix_A) - c(mpg, wt, hp)在动手分解前我们必须对这个“病人”做一次全面的体检。最关键的指标是条件数Condition Number它衡量矩阵的“病态”程度。条件数越大矩阵越接近奇异不可逆数值计算就越危险。# 计算 A^T A 的条件数这是传统最小二乘的“痛点” ATA - t(matrix_A) %*% matrix_A cond_ATA - kappa(ATA) cat(A^T A 的条件数:, cond_ATA, \n) # 输出约 1.2e04 # 计算 A 本身的条件数这是 QR 的“舒适区” cond_A - kappa(matrix_A) cat(A 本身的条件数:, cond_A, \n) # 输出约 1.1e02看到了吗A^T A的条件数比A本身高出了一百倍这正是 QR 分解要解决的核心矛盾。我们的目标就是绕过这个1.2e04的深渊直接在1.1e02的安全地带工作。4.2 手工实现 Gram-Schmidt构建 Q 矩阵现在我们抛开所有现成函数用最原始的向量运算亲手搭建 Q 矩阵。这不仅是技术练习更是对正交化思想的一次深刻内化。# 初始化 Q 矩阵与 A 同型 Q - matrix(0, nrow nrow(matrix_A), ncol ncol(matrix_A)) # 第一步处理第一列 a1 - matrix_A[, 1] q1 - a1 / sqrt(sum(a1^2)) # 单位化 Q[, 1] - q1 # 第二步处理第二列 a2 - matrix_A[, 2] # 计算 a2 在 q1 上的投影 proj_a2_q1 - sum(a2 * q1) * q1 # 减去投影得到正交分量 v2 - a2 - proj_a2_q1 # 单位化 q2 - v2 / sqrt(sum(v2^2)) Q[, 2] - q2 # 第三步处理第三列 a3 - matrix_A[, 3] # 计算 a3 在 q1 和 q2 张成的平面上的投影 proj_a3_q1 - sum(a3 * q1) * q1 proj_a3_q2 - sum(a3 * q2) * q2 v3 - a3 - proj_a3_q1 - proj_a3_q2 q3 - v3 / sqrt(sum(v3^2)) Q[, 3] - q3 # 验证 Q 的正交性Q^T %*% Q 应该是单位矩阵 QTQ - t(Q) %*% Q print(Q^T %*% Q (应为单位矩阵):) print(QTQ) # 你会看到对角线是 1非对角线是 ~1e-16证明正交性良好这段代码的每一行都在执行一个明确的几何操作。当你运行它时你看到的不再是冰冷的数字而是三个相互垂直的、长度为 1 的向量它们共同构成了一个新的坐标系。这个坐标系就是原始数据matrix_A的“骨架”。4.3 构建 R 矩阵与完整验证有了 QR 的构建就水到渠成了。根据A QR我们可以推导出R Q^T A。因为 Q 是正交的这个乘法就是将原始数据 A 投影到新的正交基 Q 上。# 构建 R 矩阵 R - t(Q) %*% matrix_A # 验证 R 是否为上三角矩阵 print(R 矩阵:) print(R) # 观察R[2,1] 和 R[3,1] 应该是极小的数~1e-16R[3,2] 同理 # 最终验证Q %*% R 是否等于原始 A reconstructed_A - Q %*% R # 计算误差 error - max(abs(matrix_A - reconstructed_A)) cat(重构误差最大值:, error, \n) # 输出约 1e-14证明分解成功此时你手中的Q和R就是一个货真价实的 QR 分解结果。虽然它的数值精度不如qr()函数但它完全是你亲手打造的、逻辑透明的“作品”。你可以清晰地看到R的第一列[R11, 0, 0]就是a1在q1方向上的分量R的第二列[R12, R22, 0]就是a2在q1和q2方向上的分量……每一个数字都有其明确的几何意义。4.4 应用实战用 QR 解决最小二乘问题现在让我们把这套“手工武器”投入到真正的战斗中——求解一个最小二乘问题。假设我们想用wt车重和hp马力来预测mpg油耗即y Xβ ε其中y是mpg列X是wt和hp组成的 32x2 矩阵。# 准备数据 y - matrix_A[, 1] # mpg 作为响应变量 X - matrix_A[, 2:3] # wt 和 hp 作为预测变量 # 对 X 进行 QR 分解注意这里 X 是 32x2所以 Q 是 32x2R 是 2x2 # 我们复用上面的 Gram-Schmidt 逻辑但只对 X 的两列操作 Q_X - matrix(0, nrow nrow(X), ncol ncol(X)) a1_X - X[, 1] q1_X - a1_X / sqrt(sum(a1_X^2)) Q_X[, 1] - q1_X a2_X - X[, 2] proj_a2_q1 - sum(a2_X * q1_X) * q1_X v2_X - a2_X - proj_a2_q1 q2_X - v2_X / sqrt(sum(v2_X^2)) Q_X[, 2] - q2_X R_X - t(Q_X) %*% X # 计算 Q^T y QTy - t(Q_X) %*% y # 现在解上三角系统 R_X %*% beta QTy # 回代法 beta - numeric(2) beta[2] - QTy[2] / R_X[2, 2] beta[1] - (QTy[1] - R_X[1, 2] * beta[2]) / R_X[1, 1] cat(QR 解得的回归系数 (wt, hp):, beta, \n) # 与 R 内置 lm() 的结果对比 model_lm - lm(y ~ X - 1) # -1 表示不加截距项与我们的 QR 设置一致 cat(lm() 解得的系数:, coef(model_lm), \n)你会发现两个结果几乎完全一致。这证明了我们手工实现的 QR 分解已经具备了实际应用的价值。它不再是教科书里的玩具而是你手中一把可以随时出鞘、解决真实问题的利剑。5. 常见问题与排错指南那些踩过的坑我都替你趟过了在无数次的调试、测试和教学中我总结了关于 QR 分解最常遇到的几类问题。它们往往不是代码语法错误而是对概念、对工具、对数据的误解。避开这些坑能帮你节省至少 80% 的无效调试时间。5.1 “我的 Q 矩阵不正交”——理解数值精度的边界这是新手最常发出的惊呼。当你计算t(Q) %*% Q期望看到一个完美的单位矩阵I结果却看到一堆1.000000e00和1.110223e-16时第一反应往往是“我哪里写错了”。注意1.110223e-16是 R 中double类型的机器精度2^-52的典型量级。它不是错误而是浮点运算的“指纹”。只要非对角线元素的绝对值小于1e-14并且对角线元素与 1 的偏差也在此量级你的 Q 矩阵就是数值上正交的。强行用round(t(Q) %*% Q, 15)去“美化”结果反而会破坏其数学一致性。排错技巧永远用max(abs(t(Q) %*% Q - diag(ncol(Q))))来量化正交性误差。如果这个值在1e-14量级恭喜你一切正常。5.2 “R 矩阵的下三角不是零”——区分算法与实现当你用 Gram-Schmidt 手写 R 矩阵时发现R[2,1]不是 0而是5.814446e-16不要慌。这恰恰证明了 Gram-Schmidt 的数值不稳定性正在起作用。而当你用qr.R(qr(X))时看到的是严格的0这是因为 LAPACK 的 Householder 实现通过精妙的算法设计将这种误差压制到了更低的量级甚至在输出时做了“四舍五入”处理。排错技巧不要拿两种不同算法的中间结果去直接比较。比较它们的最终应用效果如回归系数、重构误差才有意义。Householder 的 R 更“干净”Gram-Schmidt 的 R 更“诚实”两者各有其教学和工程价值。5.3 “为什么我的qr()结果和别人的不一样”——理解列置换PivotingR 的qr()函数有一个关键参数pivot TRUE默认。当开启列置换时qr()会自动重排X的列顺序以提升数值稳定性。这意味着qr(X)$pivot会返回一个索引向量告诉你原始列是如何被重新排列的。如果你直接用qr.Q()和qr.R()提取矩阵得到的Q和R是对应于置换后的X的。排错技巧如果你需要严格对应原始列顺序的结果务必在调用时显式指定qr(X, pivot FALSE)。否则在解读R矩阵时你会发现自己根本搞不清R[1,1]对应的是原始数据的哪个特征。5.4 “在 PySpark 里怎么用 QR”——分布式计算的思维转换在 Spark 中你无法像在单机 R 里那样直接对一个巨大的 RDD 调用qr()。分布式 QR 的核心思想是“分而治之”先将大矩阵按行切分成多个块对每个块分别进行 QR 分解得到一系列(Q_i, R_i)然后将所有R_i堆叠成一个更大的矩阵再对这个矩阵进行第二次 QR 分解最终将两次分解的Q矩阵组合起来。排错技巧不要试图在 Spark 中“复刻”单机算法。直接使用 MLlib 提供的RowMatrix类它内置了computeQR()方法。你的精力应该放在数据预处理如标准化、特征工程和结果解释上而不是底层的矩阵运算。5.5 “QR 能用来降维吗”——Q 和 R 的角色再辨析这是一个常见的误解。QR 分解本身不是一种降维技术比如 PCA 或 t-SNE。Q矩阵的列数和A的列数相同它并没有减少维度。但是R矩阵的上三角结构为我们提供了降维的“入口”。例如如果我们只保留R的前k行k n那么Q_k * R_k就是对A的一个k秩近似。这与 SVD 的U_k * Σ_k * V_k^T在功能上是等价的。排错技巧如果你想做降维QR 是一个可行的、计算成本更低的替代方案但它需要你主动“截断”R矩阵。不要指望qr()函数会自动给你一个低秩表示。6. 从 QR 到更广阔的世界它在现代机器学习栈中的位置QR 分解绝非一个孤立的、停留在教科书里的古老算法。它是现代数据科学基础设施中一根看不见的“承重柱”默默地支撑着上层无数炫目的应用。理解它在技术栈中的位置能帮你建立起更宏大的知识图谱。6.1 它是“数值计算基石”的一部分在R、Python (NumPy/SciPy)、Julia等科学计算语言中QR 分解是linalg线性代数模块的标配。它和 Cholesky 分解、LU 分解、SVD 一起构成了求解线性系统、计算矩阵特征值、进行奇异值分析的“四大金刚”。它们共同的特点是将一个复杂的、通用的矩阵问题转化为一系列结构简单、计算稳定的子问题。当你调用numpy.linalg.lstsq()时你调用的不是一个函数而是一整套经过数十年锤炼的、由 Fortran 编写的、针对不同硬件深度优化的数值计算库LAPACK/BLAS。6.2 它是“统计建模引擎”的默认驱动R 的lm()、glm()Python 的statsmodels乃至商业软件 SAS 和 Stata在进行线性回归拟合时其底层求解器的默认选项几乎无一例外是 QR 分解。原因无他**在绝大多数