数值线性方程组中的最小二乘问题解法
字数 3647 2025-12-23 10:59:52

好的，我将为您生成一个未被讲述过的重要计算数学词条。

数值线性方程组中的最小二乘问题解法

第一步：最小二乘问题的来源与定义

首先，我们从一个经典问题出发：曲线拟合。
假设我们通过实验得到了一组数据点 (x1, y1), (x2, y2), ..., (xm, ym)。我们猜测这些点近似满足一个二次多项式 y = a0 + a1*x + a2*x^2 的规律，但我们不知道系数 a0, a1, a2 是多少。

我们想找到一组“最佳”的系数，使得多项式计算出来的 y 值与实验测得的 y 值总体上最接近。如何定义“最佳”？最常用的标准是最小二乘准则：使得所有数据点上，计算值 a0 + a1*xi + a2*xi^2 与观测值 yi 之差的平方和达到最小。

将这个思想数学化。我们为每个数据点列一个方程：
a0 + a1*x1 + a2*x1^2 ≈ y1
a0 + a1*x2 + a2*x2^2 ≈ y2
...
a0 + a1*xm + a2*xm^2 ≈ ym

这可以写成一个线性方程组的形式：A a ≈ b。

  • A 是一个 m x 3 的矩阵，称为设计矩阵。第 i 行是 [1, xi, xi^2]。
  • a 是 3 x 1 的未知系数列向量 [a0; a1; a2]。
  • b 是 m x 1 的观测值列向量 [y1; y2; ...; ym]。

当 m > 3（即方程个数多于未知数个数）时，这个方程组通常是不相容的（或称为“超定的”），即不存在一个精确的解向量 a 能同时满足所有方程。此时，最小二乘解 a* 被定义为：在所有可能的向量 a 中，使得残差向量 r = b - A a 的2-范数平方（即各分量平方和）最小的那个。

用数学公式表示，最小二乘问题为：

求 a*，使得目标函数 f(a) = ||b - A a||₂² 最小化。

其中 ||·||₂ 表示向量的欧几里得范数。

第二步：几何解释与正规方程

为什么最小二乘准则在几何上很直观？我们将矩阵 A 的列向量张成一个子空间（在我们的例子中是三维空间中的一个最多三维的线性子空间，但位于 m 维空间中）。观测向量 b 是一个 m 维向量，它一般不在这个子空间里。

最小二乘解的目标，就是在子空间中寻找一个点 A a*，使得它到点 b 的欧几里得距离最短。这等价于：残差向量 r = b - A a* 必须垂直于子空间 A。

“垂直”意味着 r 与 A 的每一列都正交（点积为零）。将所有列的正交条件组合起来，就得到了一个关键的方程：

Aᵀ (b - A a*) = 0

整理后得到：

Aᵀ A a* = Aᵀ b

这个方程组被称为正规方程。Aᵀ A 是一个 n x n 的对称方阵（n 是未知数个数，此处为3），Aᵀ b 是一个 n 维向量。

核心思路的转变：求解原超定方程组 A a ≈ b 的最小二乘问题，被转化为求解一个规模更小（n x n）的正方线性方程组（正规方程）的问题。

第三步：基本数值解法——通过正规方程求解

直接求解正规方程 Aᵀ A a* = Aᵀ b 是最经典的方法。步骤如下：

  1. 计算矩阵乘法：计算 C = Aᵀ A 和 d = Aᵀ b。这需要 O(m n²) 次浮点运算。
  2. 解对称正定线性系统：求解 C a* = d。由于在大多数实际问题中（当 A 的列线性无关时），C 是对称正定矩阵，我们可以使用楚列斯基分解来高效稳定地求解：
    • 将 C 分解为 C = L Lᵀ，其中 L 是下三角矩阵。
    • 然后通过向前和向后替换求解 L y = d 和 Lᵀ a* = y。

优势：概念简单，当 n 较小且 A 条件较好时，是高效可靠的方法。

致命弱点：条件数的放大。矩阵 C = Aᵀ A 的条件数是原矩阵 A 条件数的平方，即 cond(Aᵀ A) ≈ [cond(A)]²。如果 A 本身是轻度病态的（例如，拟合多项式的基函数 1, x, x^2 在 x 值聚集时会导致 A 的列近似相关），那么 Aᵀ A 的条件数会变得非常大，使得在浮点运算中求解正规方程极不稳定，解会因舍入误差而严重失真。

第四步：更稳定的数值解法——QR分解法

为了克服正规方程法的数值不稳定性，我们需要一种不显式形成 Aᵀ A 的方法。QR分解 是解决这个问题的标准工具。

核心思想：将 m x n 矩阵 A 分解为：

A = Q R

其中：

  • Q 是一个 m x m 的正交矩阵（Qᵀ Q = I，其列向量构成 m 维空间的一组标准正交基）。
  • R 是一个 m x n 的上三角矩阵。其前 n 行组成一个上三角矩阵 R1，后 (m-n) 行全为零。

由于正交变换不改变向量的2-范数，目标函数可以改写：
||b - A a||₂² = ||Qᵀ b - Qᵀ A a||₂² = ||Qᵀ b - R a||₂²

将 Qᵀ b 也分块为 [c1; c2]，其中 c1 是前 n 行，c2 是后 m-n 行。则目标函数变为：
||b - A a||₂² = ||c1 - R1 a||₂² + ||c2||₂²

因为 ||c2||₂² 是常数，与 a 无关。所以，最小化原问题等价于最小化 ||c1 - R1 a||₂²。而当 R1 可逆时（即 A 列满秩），其最小值是0。因此，最小二乘解 a* 可以通过求解一个上三角方程组得到：

R1 a* = c1

这个过程完全避免了计算 Aᵀ A，只使用了正交变换和回代求解，其数值稳定性远高于正规方程法。

常用算法：对矩阵 A 进行 QR 分解通常使用Householder反射变换或Givens旋转来实现，它们都是数值稳定的正交变换。

第五步：处理病态问题的进阶方法——奇异值分解法

当矩阵 A 是严重病态的 或秩亏损的（即列向量线性相关或近似相关）时，QR分解法得到的解可能不稳定，甚至 R1 是奇异的。此时，最强大、最稳定的工具是奇异值分解。

核心思想：将任意 m x n 矩阵 A 分解为：

A = U Σ Vᵀ

其中：

  • U 是 m x m 的正交矩阵。
  • V 是 n x n 的正交矩阵。
  • Σ 是 m x n 的对角矩阵，其对角线元素 σ₁ ≥ σ₂ ≥ ... ≥ σᵣ ≥ 0 称为奇异值，r 是矩阵 A 的秩。

利用 SVD，原最小二乘问题的最优解（实际上是最小范数最小二乘解）具有一个简洁的表达式：

a* = V Σ⁺ Uᵀ b

其中 Σ⁺ 是一个 n x m 的对角矩阵，其对角线元素为：如果 σᵢ > 0，则对应 1/σᵢ；如果 σᵢ = 0，则对应 0。这个矩阵被称为 Σ 的伪逆。

优势：

  1. 稳定性：SVD 是数值线性代数中最稳定的分解之一，能精确揭示矩阵的秩和病态程度。
  2. 处理秩亏损：当 σᵢ 非常小（接近机器精度）时，我们可以将其视为零，从而截断 SVD。在解的表达式中，将对应的 1/σᵢ 置为0，这相当于用低秩矩阵近似了原矩阵 A，过滤掉了由小奇异值引起的数值噪声，从而得到一个正则化的、数值稳定的解。这种技术称为截断奇异值分解正则化。

代价：SVD 的计算成本是三种方法中最高的，约为 O(m n²)，但它提供了最深刻的数值洞察和最稳定的解。

总结

求解最小二乘问题 min ||b - A a||₂ 的数值方法是一个从简单到稳健、从常规到处理病态的完整谱系：

  1. 正规方程法 (Aᵀ A a = Aᵀ b)：思想基础，适合小规模、良条件问题。不推荐作为通用方法，因其数值不稳定性。
  2. QR分解法 (A = Q R, 解 R1 a = c1)：标准方法。通过正交变换避免了条件数放大，对绝大多数问题稳定可靠。
  3. 奇异值分解法 (A = U Σ Vᵀ, a* = V Σ⁺ Uᵀ b)：终极武器。计算成本最高，但能完美处理秩亏损和严重病态问题，并提供正则化解。是分析问题数值特性、进行模型诊断的黄金标准。

在实际的数值计算软件（如 MATLAB、NumPy/SciPy）中，当你调用 A \ b 或类似的最小二乘求解器时，库函数通常会根据矩阵 A 的性质自动在 QR 分解和 SVD 等方法之间选择最合适、最高效的算法。

好的，我将为您生成一个未被讲述过的重要计算数学词条。 数值线性方程组中的最小二乘问题解法 第一步：最小二乘问题的来源与定义 首先，我们从一个经典问题出发： 曲线拟合 。 假设我们通过实验得到了一组数据点 (x1, y1), (x2, y2), ..., (xm, ym) 。我们猜测这些点近似满足一个二次多项式 y = a0 + a1*x + a2*x^2 的规律，但我们不知道系数 a0, a1, a2 是多少。 我们想找到一组“最佳”的系数，使得多项式计算出来的 y 值与实验测得的 y 值总体上最接近。如何定义“最佳”？最常用的标准是 最小二乘准则 ：使得所有数据点上，计算值 a0 + a1*xi + a2*xi^2 与观测值 yi 之差的 平方和 达到最小。 将这个思想数学化。我们为每个数据点列一个方程： a0 + a1*x1 + a2*x1^2 ≈ y1 a0 + a1*x2 + a2*x2^2 ≈ y2 ... a0 + a1*xm + a2*xm^2 ≈ ym 这可以写成一个 线性方程组 的形式： A a ≈ b 。 A 是一个 m x 3 的矩阵，称为 设计矩阵 。第 i 行是 [1, xi, xi^2] 。 a 是 3 x 1 的未知系数列向量 [a0; a1; a2] 。 b 是 m x 1 的观测值列向量 [y1; y2; ...; ym] 。 当 m > 3 （即方程个数多于未知数个数）时，这个方程组通常是 不相容的 （或称为“超定的”），即不存在一个精确的解向量 a 能同时满足所有方程。此时，最小二乘解 a * 被定义为：在所有可能的向量 a 中，使得 残差向量 r = b - A a 的 2-范数平方 （即各分量平方和）最小的那个。 用数学公式表示，最小二乘问题为： 求 a * ，使得目标函数 f(a) = ||b - A a||₂² 最小化。 其中 ||·||₂ 表示向量的欧几里得范数。 第二步：几何解释与正规方程 为什么最小二乘准则在几何上很直观？我们将矩阵 A 的列向量张成一个 子空间 （在我们的例子中是三维空间中的一个最多三维的线性子空间，但位于 m 维空间中）。观测向量 b 是一个 m 维向量，它一般不在这个子空间里。 最小二乘解的目标，就是在子空间中寻找一个点 A a * ，使得它到点 b 的 欧几里得距离最短 。这等价于：残差向量 r = b - A a * 必须垂直于子空间 A 。 “垂直”意味着 r 与 A 的 每一列 都正交（点积为零）。将所有列的正交条件组合起来，就得到了一个关键的方程： Aᵀ (b - A a* ) = 0 整理后得到： Aᵀ A a* = Aᵀ b 这个方程组被称为 正规方程 。 Aᵀ A 是一个 n x n 的对称方阵（ n 是未知数个数，此处为3）， Aᵀ b 是一个 n 维向量。 核心思路的转变 ：求解原超定方程组 A a ≈ b 的最小二乘问题，被转化为求解一个规模更小（ n x n ）的 正方线性方程组 （正规方程）的问题。 第三步：基本数值解法——通过正规方程求解 直接求解正规方程 Aᵀ A a* = Aᵀ b 是最经典的方法。步骤如下： 计算矩阵乘法 ：计算 C = Aᵀ A 和 d = Aᵀ b 。这需要 O(m n²) 次浮点运算。 解对称正定线性系统 ：求解 C a* = d 。由于在大多数实际问题中（当 A 的列线性无关时）， C 是对称正定矩阵，我们可以使用 楚列斯基分解 来高效稳定地求解： 将 C 分解为 C = L Lᵀ ，其中 L 是下三角矩阵。 然后通过向前和向后替换求解 L y = d 和 Lᵀ a* = y 。 优势 ：概念简单，当 n 较小且 A 条件较好时，是高效可靠的方法。 致命弱点 ： 条件数 的放大。矩阵 C = Aᵀ A 的条件数是原矩阵 A 条件数的平方，即 cond(Aᵀ A) ≈ [cond(A)]² 。如果 A 本身是轻度病态的（例如，拟合多项式的基函数 1, x, x^2 在 x 值聚集时会导致 A 的列近似相关），那么 Aᵀ A 的条件数会变得非常大，使得在浮点运算中求解正规方程 极不稳定 ，解会因舍入误差而严重失真。 第四步：更稳定的数值解法——QR分解法 为了克服正规方程法的数值不稳定性，我们需要一种不显式形成 Aᵀ A 的方法。 QR分解 是解决这个问题的标准工具。 核心思想 ：将 m x n 矩阵 A 分解为： A = Q R 其中： Q 是一个 m x m 的 正交矩阵 （ Qᵀ Q = I ，其列向量构成 m 维空间的一组标准正交基）。 R 是一个 m x n 的 上三角矩阵 。其前 n 行组成一个上三角矩阵 R1 ，后 (m-n) 行全为零。 由于正交变换不改变向量的2-范数，目标函数可以改写： ||b - A a||₂² = ||Qᵀ b - Qᵀ A a||₂² = ||Qᵀ b - R a||₂² 将 Qᵀ b 也分块为 [c1; c2] ，其中 c1 是前 n 行， c2 是后 m-n 行。则目标函数变为： ||b - A a||₂² = ||c1 - R1 a||₂² + ||c2||₂² 因为 ||c2||₂² 是常数，与 a 无关。所以，最小化原问题等价于 最小化 ||c1 - R1 a||₂² 。而当 R1 可逆时（即 A 列满秩），其最小值是0。因此，最小二乘解 a * 可以通过求解一个上三角方程组得到： R1 a* = c1 这个过程完全避免了计算 Aᵀ A ，只使用了正交变换和回代求解，其数值稳定性远高于正规方程法。 常用算法 ：对矩阵 A 进行 QR 分解通常使用 Householder反射变换 或 Givens旋转 来实现，它们都是数值稳定的正交变换。 第五步：处理病态问题的进阶方法——奇异值分解法 当矩阵 A 是 严重病态的 或 秩亏损的 （即列向量线性相关或近似相关）时，QR分解法得到的解可能不稳定，甚至 R1 是奇异的。此时，最强大、最稳定的工具是 奇异值分解 。 核心思想 ：将任意 m x n 矩阵 A 分解为： A = U Σ Vᵀ 其中： U 是 m x m 的正交矩阵。 V 是 n x n 的正交矩阵。 Σ 是 m x n 的对角矩阵，其对角线元素 σ₁ ≥ σ₂ ≥ ... ≥ σᵣ ≥ 0 称为奇异值， r 是矩阵 A 的秩。 利用 SVD，原最小二乘问题的最优解（实际上是最小范数最小二乘解）具有一个简洁的表达式： a* = V Σ⁺ Uᵀ b 其中 Σ⁺ 是一个 n x m 的对角矩阵，其对角线元素为：如果 σᵢ > 0 ，则对应 1/σᵢ ；如果 σᵢ = 0 ，则对应 0 。这个矩阵被称为 Σ 的伪逆。 优势 ： 稳定性 ：SVD 是数值线性代数中最稳定的分解之一，能精确揭示矩阵的秩和病态程度。 处理秩亏损 ：当 σᵢ 非常小（接近机器精度）时，我们可以将其视为零，从而 截断 SVD。在解的表达式中，将对应的 1/σᵢ 置为0，这相当于用低秩矩阵近似了原矩阵 A ，过滤掉了由小奇异值引起的数值噪声，从而得到一个 正则化的、数值稳定的解 。这种技术称为 截断奇异值分解正则化 。 代价 ：SVD 的计算成本是三种方法中最高的，约为 O(m n²) ，但它提供了最深刻的数值洞察和最稳定的解。 总结 求解最小二乘问题 min ||b - A a||₂ 的数值方法是一个从简单到稳健、从常规到处理病态的完整谱系： 正规方程法 ( Aᵀ A a = Aᵀ b ) ：思想基础，适合小规模、良条件问题。 不推荐 作为通用方法，因其数值不稳定性。 QR分解法 ( A = Q R , 解 R1 a = c1 ) ： 标准方法 。通过正交变换避免了条件数放大，对绝大多数问题稳定可靠。 奇异值分解法 ( A = U Σ Vᵀ , a* = V Σ⁺ Uᵀ b ) ： 终极武器 。计算成本最高，但能完美处理秩亏损和严重病态问题，并提供正则化解。是分析问题数值特性、进行模型诊断的黄金标准。 在实际的数值计算软件（如 MATLAB、NumPy/SciPy）中，当你调用 A \ b 或类似的最小二乘求解器时，库函数通常会根据矩阵 A 的性质自动在 QR 分解和 SVD 等方法之间选择最合适、最高效的算法。