最小二乘法的正规方程推导及条件数对解稳定性的影响

构建最小二乘问题的基本模型。我们有一个线性方程组 $Ax=b$，其中 $A$ 是一个 $m \times n$ 的矩阵（通常 $m > n$，即方程个数多于未知数个数），$b$ 是一个 $m$ 维的观测向量。由于数据可能存在噪声或模型不完全精确，这个方程组很可能无精确解。最小二乘法的目标是寻找一个 $n$ 维向量 $x$，使得残差向量 $r = Ax - b$ 的长度（即 $\|r\|_2^2$）最小。

定义目标函数。将残差的平方和作为优化的目标，即：
$$ J(x) = \|Ax - b\|_2^2 = (Ax-b)^T(Ax-b) $$
我们的任务就是最小化这个二次函数 $J(x)$。

1. 推导正规方程

展开目标函数 $J(x)$。利用向量和矩阵的转置运算规则：
$$ J(x) = (Ax-b)^T(Ax-b) = (x^TA^T - b^T)(Ax - b) = x^TA^TAx - x^TA^Tb - b^TAx + b^Tb $$
注意到 $x^TA^Tb$ 是一个标量（1x1矩阵），而标量的转置等于其自身。因为 $(x^TA^Tb)^T = b^TAx$，所以 $x^TA^Tb = b^TAx$。因此，目标函数可简化为：
$$ J(x) = x^TA^TAx - 2b^TAx + b^Tb $$

对目标函数求导。利用矩阵微积分中关于向量求导的公式：$\frac{\partial (x^TAx)}{\partial x} = 2Ax$ (当 $A$ 对称时)，$\frac{\partial (b^TAx)}{\partial x} = A^Tb$。这里 $A^TA$ 是对称矩阵。因此，对 $J(x)$ 关于 $x$ 求梯度：
$$ \nabla J(x) = \frac{\partial J(x)}{\partial x} = 2A^TAx - 2A^Tb $$

设置导数为零以找到极小值点。要使 $J(x)$ 最小，其梯度应为零向量，即：
$$ 2A^TAx - 2A^Tb = 0 $$
两边同时除以 2 并移项，得到最终的 正规方程：
$$ A^TAx = A^Tb $$

求解正规方程。如果矩阵 $A^TA$ 是可逆的（即满秩），则最小二乘解的解析形式为：
$$ x = (A^TA)^{-1}A^Tb $$

2. 理解条件数及其计算

认识矩阵的条件数。条件数是衡量一个矩阵（或线性方程组）对输入数据的小扰动有多敏感的数值指标。对于一个可逆矩阵 $M$，其条件数 $\kappa(M)$ 定义为：
$$ \kappa(M) = \|M\| \cdot \|M^{-1}\| $$
其中 $\|\cdot\|$ 通常取矩阵的 2-范数（即最大奇异值）。条件数 $\kappa(M)$ 总是大于等于 1。

计算 $A^TA$ 的条件数。在最小二乘问题中，我们关心的是求解正规方程时系统矩阵 $A^TA$ 的条件数，记作 $\kappa(A^TA)$。根据条件数的性质，有 $\kappa(A^TA) = \kappa(A)^2$。这意味着如果原始数据矩阵 $A$ 的条件数已经很大，那么 $A^TA$ 的条件数会急剧增大，变得非常“病态”。

实际计算步骤（以计算 $\kappa(A)$ 为例）：

1. 输入矩阵 $A$。
2. 计算矩阵 $A$ 的奇异值分解 (SVD)：$A = U\Sigma V^T$，其中 $\Sigma$ 是一个对角阵，其对角线元素 $\sigma_1 \geq \sigma_2 \geq ... \geq \sigma_n > 0$ 是 $A$ 的奇异值。
3. 获取最大奇异值 $\sigma_{\text{max}} = \sigma_1$ 和最小奇异值 $\sigma_{\text{min}} = \sigma_n$。
4. 计算矩阵 $A$ 的 2-范数条件数：$\kappa(A) = \sigma_{\text{max}} / \sigma_{\text{min}}$。
5. 得到 $A^TA$ 的条件数：$\kappa(A^TA) = (\kappa(A))^2$。

3. 条件数对解稳定性的影响

判断数值稳定性。条件数 $\kappa(A^TA)$ 的大小直接决定了用正规方程求解最小二乘问题的数值稳定性。

• 当 $\kappa(A^TA)$ 较小（接近 1）时，矩阵 $A^TA$ 是“良态”的。此时，观测向量 $b$ 中的微小误差（噪声）对解 $x$ 的影响也会很小。求解过程稳定可靠。

• 当 $\kappa(A^TA)$ 非常大（例如 $10^{10}$ 或更大）时，矩阵 $A^TA$ 是“病态”的。此时，$b$ 中即使只有舍入误差级别的微小扰动，也会在求解 $x = (A^TA)^{-1}A^Tb$ 的过程中被剧烈放大，导致计算出的 $x$ 与真实解相差甚远，结果完全不可信。

具体表现为：

1. 求解困难：在用高斯消元法或 Cholesky 分解求解 $A^TAx = A^Tb$ 时，主元可能非常小，导致数值溢出或极大的舍入误差。
2. 解剧烈变化：对 $b$ 做微小改动（例如改变 $10^{-6}$），解 $x$ 的某些分量可能变化巨大。

4. 提高数值稳定性的实践指南

当检测到问题病态（条件数很大）时，避免直接使用正规方程。采用以下数值更稳定的方法：

1. 使用 QR 分解。对矩阵 $A$ 进行 QR 分解：$A = QR$，其中 $Q$ 是正交矩阵，$R$ 是上三角矩阵。最小二乘问题等价于求解 $Rx = Q^Tb$。由于 $Q$ 是正交矩阵，它不改变向量的 2-范数，因此这个变换是数值稳定的。求解 $Rx = Q^Tb$ 仅需简单的回代，避免了形成病态矩阵 $A^TA$。

2. 使用奇异值分解 (SVD)。计算 $A$ 的 SVD：$A = U\Sigma V^T$。最小二乘解由 $x = V\Sigma^{\dagger}U^Tb$ 给出，其中 $\Sigma^{\dagger}$ 是 $\Sigma$ 的伪逆（将非零奇异值取倒数）。SVD 是处理病态问题最稳健的方法，它能清晰地揭示问题的数值秩，并允许通过截断小奇异值进行正则化。

3. 特征值/奇异值检查。在生成 $A$ 和 $b$ 后，计算 $\kappa(A)$。如果发现 $\kappa(A) > 10^6$，就需要警惕，并优先考虑使用 QR 或 SVD 方法。

4. 数据中心化与缩放。有时，对数据矩阵 $A$ 的列进行中心化（减去均值）和缩放（除以标准差），可以显著降低条件数 $\kappa(A)$，从而改善正规方程的数值性质。

评估你的解。无论使用哪种方法，都应计算残差 $r = Ax - b$ 并检查其范数是否合理。对于病态问题，可以尝试对解 $x$ 施加微小扰动，观察残差的变化，以感受解的敏感性。