误差分析与超定方程组¶

参考：计算方法 (00330050) 第 4 讲《解线性方程组的直接法 (误差分析, 超定方程组)》

目录¶

• 简介
• 向量范数与矩阵范数
• 误差分析
• 条件数的估计
• 条件数的再思考
• 超定方程组与最小二乘

简介：扰动从哪里来¶

在计算机中，实数以浮点数形式存储，因此会产生舍入误差。
对数学上精确给定的线性方程组

输入计算机后通常变为

其中：

• \(\delta A\) 是系数矩阵扰动
• \(\delta b\) 是右端项扰动
• \(\delta x\) 是解的变化（误差传播结果）

如果很小的输入扰动会导致很大的解变化，则称该问题为 病态问题（矩阵是 病态矩阵）；反之称为 良态问题。

例：病态方程组¶

考虑

原解为

若右端项微扰为 \(\delta b = [0,\, 0.0001]^T\)，则新解变成

所以

说明右端极小变化引起解的大幅变化，该系统病态。

向量范数与矩阵范数¶

向量范数定义¶

向量范数 \(\|x\|\) 需要满足：

• 正定性：\(\|x\| \ge 0\)，且 \(\|x\|=0 \Leftrightarrow x=0\)
• 齐次性：\(\|\alpha x\| = |\alpha| \|x\|\)
• 三角不等式：\(\|x+y\| \le \|x\|+\|y\|\)

常见向量范数：

范数等价¶

在有限维空间 \(\mathbb{R}^n\) 中，任意两个范数等价。典型不等式：

矩阵范数定义¶

矩阵范数 \(\|A\|\) 除了正定性、齐次性、三角不等式外，还需满足相容性：

常见矩阵范数：

诱导矩阵范数与谱半径¶

由向量范数诱导的矩阵范数定义为

谱半径定义：

这里 \(\lambda_i\) 为矩阵的特征值

并有一般结论

如果对 \(\forall A \in \mathbb{R}^{n \times n}\) 与 \(\forall x \in \mathbb{R}^n\)，都有

那么，称式中的矩阵范数与向量范数 相容。

误差分析：扰动如何放大¶

仅右端项扰动¶

若

则

从而

结合 \(b=Ax\) 可得

最终得到相对误差估计

即解 \(x\) 的相对误差不会超过 \(b\) 的相对误差的 \(\|A\|\|A^{-1}\|\) 倍

仅系数矩阵扰动¶

若

整理得

进一步可得

当 \(\|A^{-1}\|\,\|\delta A\|<1\) 时成立。

综合扰动结论（先验估计）¶

若同时有 \(\delta A,\delta b\)，并且 \(\|A^{-1}\|\,\|\delta A\|<1\)，则

条件数¶

定义（对非奇异矩阵）：

按不同范数记为

条件数性质¶

• \(\mathrm{Cond}(A)\ge 1\)
• 若 \(A\) 为正交矩阵（2 范数下），则 \(\mathrm{Cond}_2(A)=1\)，也就是正交矩阵是一种典型的良态矩阵。
• 对任意 \(\alpha \ne 0\)，有 \(\mathrm{Cond}(\alpha A)=\mathrm{Cond}(A)\)
• 在 2 范数下：

不同范数下的条件数在有限维空间中等价。

条件数的估计与再思考¶

如何估计条件数¶

直接算 \(A^{-1}\) 代价高。下面给出基于 LU分解 求逆的算法及复杂度分析。

逆矩阵计算算法¶

对于非奇异矩阵 \(A \in \mathbb{R}^{n \times n}\)，利用 LU 分解 \(A = LU\) 求逆：

• 先对 \(A\) 进行 LU 分解（带部分主元），得到 \(PA = LU\)
• 然后逐列求解 \(AX = I\)，即求解 \(n\) 个线性方程组

其中 \(e_i\) 是第 \(i\) 个标准基向量。

算法复杂度分析¶

因此直接求逆的总复杂度约为 \(\frac{8}{3}n^3\)，远高于求解单个方程组的 \(\frac{2}{3}n^3\)。

算法伪代码¶

输入: A ∈ R^(n×n) (非奇异)
输出: A^(-1)

1. 对 A 进行 LU 分解（带部分主元）:  PA = LU
2. 初始化 X = 0 ∈ R^(n×n)
3. for j = 1 to n:
4.     令 e_j = [0,...,1,...,0]^T  (第j个为1)
5.     解 Ly = P*e_j      (前代)
6.     解 Ux_j = y       (回代)
7.     X(:,j) = x_j
8. return X

import numpy as np
from scipy.linalg import lu_factor, lu_solve

def matrix_inverse_lu(A):
    """
    利用 LU 分解求逆矩阵

    参数:
        A: n×n 非奇异矩阵
    返回:
        A^{-1}
    """
    n = A.shape[0]
    # LU 分解（带部分主元）
    lu, piv = lu_factor(A)

    # 构造单位矩阵的各列
    X = np.zeros((n, n))
    for j in range(n):
        e_j = np.zeros(n)
        e_j[j] = 1.0
        # 解 Ax = e_j
        X[:, j] = lu_solve((lu, piv), e_j)

    return X

# 示例
if __name__ == "__main__":
    A = np.array([[4, 3],
                  [6, 3]], dtype=float)
    A_inv = matrix_inverse_lu(A)
    print("A^{-1} =")
    print(A_inv)
    print("验证: A @ A^{-1} =")
    print(A @ A_inv)

#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <Eigen/LU>

using namespace Eigen;
using namespace std;

MatrixXd matrixInverseLU(const MatrixXd& A) {
    """
    利用 LU 分解求逆矩阵

    参数:
        A: n×n 非奇异矩阵
    返回:
        A^{-1}
    """
    int n = A.rows();

    // LU 分解（带部分主元）
    PartialPivLU<MatrixXd> lu(A);

    // 求解 n 个右端项（单位矩阵各列）
    MatrixXd I = MatrixXd::Identity(n, n);
    MatrixXd A_inv = lu.solve(I);

    return A_inv;
}

int main() {
    MatrixXd A(2, 2);
    A << 4, 3,
         6, 3;

    MatrixXd A_inv = matrixInverseLU(A);

    cout << "A^{-1} =" << endl << A_inv << endl;
    cout << "验证: A * A^{-1} =" << endl << A * A_inv << endl;

    return 0;
}

条件数估计的替代方法¶

由

可得

通过选取合适 \(d\)，让 \(\|y\|\) 尽量大，可以构造 \(\|A^{-1}\|\) 的下界估计。实际计算中常用 Hager 估计法 或 条件数迭代估计，避免显式求逆。

条件数的几何意义¶

设 \(S\) 为奇异矩阵集合。对诱导范数有：

这说明：

• 条件数大 \(\Rightarrow\) \(A\) 离奇异矩阵集合很近（病态）
• 条件数小 \(\Rightarrow\) \(A\) 离奇异集合较远（良态）

超定方程组与最小二乘¶

当方程数 \(m\) 大于未知数数 \(n\) 时，得到线性超定方程组：

一般情况下无精确解，因此转化为优化问题。

残差与目标函数¶

定义残差（余值）

残差向量（余向量）

由于一个向量本身不具备大小的可比性，因此，我们需要引入一个实函数 \(F(x)\)，使得

定义最小二乘目标函数

正规方程（Normal Equation）¶

为求 \(F(x) = (b-Ax)^T(b-Ax)\) 的极小值，令梯度 \(\nabla F(x)=0\)。

详细推导过程：

第一步：展开目标函数

令残差向量 \(r = b - Ax\)，则：

第二步：展开矩阵乘法

注释 1：注意到 \(b^T A x\) 是一个标量（\(1 \times 1\) 矩阵），而标量的转置等于自身，即 \((b^T A x)^T = x^T A^T b = b^T A x\)。因此中间两项相等：\(b^T A x = x^T A^T b\)。

第三步：合并并简化

第四步：对 \(x\) 求梯度

利用矩阵求导的基本公式：

• \(\frac{\partial}{\partial x}(c^T x) = c\)（线性项）
• \(\frac{\partial}{\partial x}(x^T Q x) = 2Qx\)（二次项，当 \(Q\) 对称时）

注释 2：这里 \(A^T A\) 是对称矩阵（因为 \((A^T A)^T = A^T A\)），所以二次项求导直接得到 \(2A^T A x\)。

逐项求导：

第五步：令梯度为零

约去因子 2，得到 正规方程：

第六步：解的存在性与唯一性

若 \(\mathrm{rank}(A)=n\)（\(A\) 列满秩），则：

• \(A^T A\) 是 \(n \times n\) 对称矩阵
• \(A^T A\) 正定（因为对任意非零 \(y\)，\(y^T A^T A y = \|Ay\|_2^2 > 0\)）
• 故 \(A^T A\) 可逆，解唯一：\(x = (A^T A)^{-1} A^T b\)

小结¶

• 误差分析核心：看输入扰动如何通过 \(A^{-1}\) 放大到解误差
• 条件数是衡量问题敏感性的核心量
• 条件数可解释为“到奇异矩阵集合的相对距离倒数”
• 超定问题通常用最小二乘，最终化为正规方程 \(A^T A x=A^T b\)