计算机科学的数学基础-03：线性系统与矩阵分解

线性系统

线性系统写作：$A\vec{x}=\vec{b},\quad A\in\mathbb{R}^{m\times n},\quad\vec{x}\in\mathbb{R}^{n},\quad\vec{b}\in\mathbb{R}^{m}$

每一行可以理解为一个线性约束。求解线性系统，就是寻找同时满足所有行约束的向量 $\vec{x}$。

特殊矩阵

对角矩阵最容易求解：

$$\begin{bmatrix} d_1 & 0 & \cdots & 0\\ 0 & d_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & d_n \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix} = \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n \end{bmatrix}$$

只要 $d_i\ne 0$，就有 $x_i=\frac{b_i}{d_i}$。

正交矩阵满足：

$$A^{-1}=A^T$$

因此：

$$\vec{x}=A^{-1}\vec{b}=A^T\vec{b}$$

上三角矩阵 $U\vec{x}=\vec{b}$ 可以从最后一行开始后向代入；下三角矩阵 $L\vec{x}=\vec{b}$ 可以从第一行开始前向代入。正交矩阵在数值计算中重要，因为转置比一般求逆更便宜，也通常更稳定。

宽 / 高矩阵

宽矩阵指 $m<n$，方程数量少于未知量数量。若系统可解，通常有自由变量，因此不可能有唯一解。

例如：

$$\begin{bmatrix} 2 & 3 & -1\\ 1 & -1 & 4 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix} 5\\ 11 \end{bmatrix}$$

它至少有两个解：

$$\begin{bmatrix} 1\\ 2\\ 3 \end{bmatrix}, \qquad \begin{bmatrix} 7.6\\ -3.4\\ 0 \end{bmatrix}$$

解集应写成“一个特解加零空间方向”的形式：

$$\vec{x} = \begin{bmatrix} 1\\ 2\\ 3 \end{bmatrix} +t \left( \begin{bmatrix} 7.6\\ -3.4\\ 0 \end{bmatrix} - \begin{bmatrix} 1\\ 2\\ 3 \end{bmatrix} \right)$$

高矩阵指 $m>n$，方程数量多于未知量数量。它可能有解，也可能因为约束矛盾而无解。更一般地，对任意高矩阵 $A$，都存在某些右端项 $\vec{b}_0$，使得 $A\vec{x}=\vec{b}_0$ 不可解。

默认关注方阵 $A\in\mathbb{R}^{n\times n}$，并且先假设 $A$ 可逆。此时对任意 $\vec{b}$，系统都有唯一解。

矩阵向量乘

矩阵向量乘法返回 $\vec{b}=A\vec{x}$。按行计算时：

Multiply(A, x):
    # A 有 m 行 n 列，x 有 n 个分量。
    b <- 0

    # 第 i 个输出分量是第 i 行和 x 的内积。
    for i <- 1, 2, ..., m:
        for j <- 1, 2, ..., n:
            b_i <- b_i + a_ij * x_j

    return b

也可以按列计算，把 $x_j$ 看作第 $j$ 列的系数：

Multiply(A, x):
    # 逐列累加每一列对 b 的贡献。
    b <- 0

    for j <- 1, 2, ..., n:
        for i <- 1, 2, ..., m:
            b_i <- b_i + a_ij * x_j

    return b

两种写法数学上等价，实际性能取决于矩阵的存储方式和缓存访问模式。

数值表示与误差

二进制表示

整数部分：

$$(b_\ell b_{\ell-1}\cdots b_0)_2=\sum_{i=0}^{\ell}b_i2^i,\qquad b_i\in\{0,1\}$$

小数部分：

$$(0.b_1b_2\cdots b_k)_2=\sum_{i=1}^{k}b_i2^{-i}$$

很多十进制小数不能用有限二进制位精确表示，例如 $\frac{1}{3}=(0.0101010101\cdots)_2$。因此，浮点计算必须处理舍入误差。

浮点数与容差

浮点数不应直接用等号比较。数学上相等的表达式，经过有限精度运算后可能只是在误差范围内接近。

不推荐：

double x = 1.0;
double y = x / 3.0;

// y * 3.0 不一定能精确回到 x。
if (x == y * 3.0)
    cout << "They_are_equal!";
else
    cout << "They_are_NOT_equal!";

推荐使用容差：

double x = 1.0;
double y = x / 3.0;
const double epsilon = numeric_limits<double>::epsilon();

// 判断两个浮点结果是否足够接近。
if (fabs(x - y * 3.0) < epsilon)
    cout << "They_are_equal!";
else
    cout << "They_are_NOT_equal!";

容差大小要结合问题尺度选择。过小会把数值等价误判为不等；过大会掩盖真实误差。

定点数

定点数固定小数点位置。若有 $k$ 位小数和 $\ell+1$ 位整数部分，总位数为 $k+\ell+1$。

位权结构为：

$$\begin{array}{cccccccc} 1 & 1 & \cdots & 0 & 0 & \cdots & 1 & 1\\ 2^\ell & 2^{\ell-1} & \cdots & 2^0 & 2^{-1} & \cdots & 2^{-k+1} & 2^{-k} \end{array}$$

定点加法可以复用整数加法：$a+b=(a\cdot 2^k+b\cdot 2^k)\cdot 2^{-k}$。

定点数实现简单，但动态范围有限。乘除法会改变数量级，例如 $(0.1)_2\times(0.1)_2=(0.01)_2$。如果只保留一位小数，结果会被近似为 $(0.0)_2$。

浮点数格式

科学记数法存储有效数字和指数：

$$\pm\left(d_0+d_1b^{-1}+d_2b^{-2}+\cdots+d_{p-1}b^{1-p}\right)\times b^e$$

其中：

• $b$ 是基数。
• $p$ 是精度，即有效数字位数。
• $e\in[L,U]$ 是指数范围。

单精度浮点数共有 32 位：

• 符号位 $S$：1 位。
• 偏置指数 $E$：8 位，偏置量为 127。
• 尾数 $F$：23 位，规格化数省略最高位 1。

若 $c=0.2998\times 10^9$，则 $c=(0001\,0001\,1101\,1110\,1001\,0101\,1100\,0000)_2$。

写成规格化二进制：

$$c=1.001\,1101\,1110\,1001\,0101\,1100\,0000\times 2^{28}$$

因此 $S=0$，$E=28+127=155_{10}=(10011011)_2$。尾数 $F$ 是去掉最高位 1 后的 23 位。

前向 / 后向误差与条件数

前向误差衡量近似解和真实解的差异：

$$\text{forward error}=|\hat{x}-x|$$

后向误差衡量需要把原问题改变多少，才能使当前近似解成为扰动后问题的精确解。

以根查找为例，真实根 $x_0$ 满足 $f(x_0)=0$。若算法输出 $x_{\mathrm{est}}$，则前向误差是 $|x_{\mathrm{est}}-x_0|$。但 $x_0$ 通常未知，无法直接计算。可计算的替代量是残差 $|f(x_{\mathrm{est}})-f(x_0)|=|f(x_{\mathrm{est}})|$。

条件数描述前向误差和后向误差之间的放大关系：

$$\text{condition number}=\frac{\text{forward error}}{\text{backward error}}$$

若小后向误差必然导致小前向误差，问题是良态的；否则问题是病态的。

对线性方程 $ax=b,\ x_0=\frac{b}{a}$，条件数为：

$$\frac{|x-x_0|}{|a(x-x_0)|}=\frac{1}{|a|}$$

当 $|a|$ 很小时，问题对扰动敏感。

对根查找问题，若近似解为 $x+\varepsilon$：

$$\frac{\text{forward error}}{\text{backward error}}=\frac{|(x+\varepsilon)-x|}{|f(x+\varepsilon)-f(x)|}\approx\frac{1}{|f'(x)|}$$

根附近 $|f'(x)|$ 越小，同样的残差可能对应越大的根位置误差。

数量级问题

数值软件中常见风险包括：

• 除以非常小的数。
• 两个非常接近的数相减。
• 中间结果溢出或下溢。

以二范数为例，直接累加平方可能溢出：

double normSquared = 0;

// x[i] 很大时，x[i] * x[i] 可能先溢出。
for (int i = 0; i < n; i++)
    normSquared += x[i] * x[i];

return sqrt(normSquared);

更稳妥的做法是先缩放：

double maxElement = epsilon;

// 找到最大绝对值作为缩放因子。
for (int i = 0; i < n; i++)
    maxElement = max(maxElement, fabs(x[i]));

double normSquared = 0;

// 缩放后再平方求和。
for (int i = 0; i < n; i++) {
    double scaled = x[i] / maxElement;
    normSquared += scaled * scaled;
}

return sqrt(normSquared) * maxElement;

对应公式为：

$$\|\vec{x}\|_2=\sqrt{\sum_i x_i^2}=\max_i|x_i|\sqrt{\sum_i\left(\frac{x_i}{\max_j|x_j|}\right)^2}$$

行操作的矩阵表示

高斯消元使用三类基本行操作：

1. 交换行。
2. 缩放行。
3. 把一行的倍数加到另一行。

这些操作都可以写成左乘矩阵。

置换矩阵

置换是一个双射 $\sigma:\{1,\dots,m\}\to\{1,\dots,m\}$。

置换矩阵为：

$$P_\sigma= \begin{bmatrix} - & \vec{e}_{\sigma(1)}^T & -\\ - & \vec{e}_{\sigma(2)}^T & -\\ & \vdots &\\ - & \vec{e}_{\sigma(m)}^T & - \end{bmatrix}$$

左乘 $P_\sigma$ 会重排矩阵的行。

行缩放矩阵

把第 $k$ 行乘以 $a_k$，可以左乘对角矩阵：

$$S_a= \begin{bmatrix} a_1 & 0 & \cdots & 0\\ 0 & a_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & a_m \end{bmatrix}$$

若每个 $a_k\ne 0$，则 $S_a^{-1}=S_{1/a}$。

消元矩阵

矩阵 $\vec{e}_\ell\vec{e}_k^T$ 只有第 $\ell$ 行第 $k$ 列为 1，其余位置为 0。左乘时，$\vec{e}_\ell\vec{e}_k^TA$ 会把 $A$ 的第 $k$ 行复制到第 $\ell$ 行，其余行为 0。

若要把第 $k$ 行的 $c$ 倍加到第 $\ell$ 行，可以左乘：

$$M=I_{n\times n}+c\vec{e}_\ell\vec{e}_k^T$$

当 $k\ne \ell$ 时，$\vec{e}_k^T\vec{e}_\ell=0$，所以：

$$M^{-1}=I_{n\times n}-c\vec{e}_\ell\vec{e}_k^T$$

高斯消元

高斯消元把增广矩阵 $(A\mid \vec{b})$ 变为上三角系统 $U\vec{x}=\vec{y}$，然后通过后向代入求解。

前向消元

前向消元的核心步骤：

1. 选择主元。
2. 必要时进行行置换。
3. 用主元行消去主元下方元素。
4. 在剩余子矩阵中重复。

伪代码：

ForwardElimination(A, b):
    # 复制输入，避免覆盖原矩阵。
    U, y <- A, b

    for p <- 1, 2, ..., n:
        optionally pivot

        # 将主元归一化。
        s <- 1 / u_pp
        y_p <- s * y_p
        for c <- p, ..., n:
            u_pc <- s * u_pc

        # 消去主元下方的元素。
        for r <- p + 1, ..., n:
            s <- -u_rp
            y_r <- y_r + s * y_p
            for c <- p, ..., n:
                u_rc <- u_rc + s * u_pc

    return U, y

若主元为 0 或过小，必须选主元，否则算法可能失败或严重不稳定。

后向代入

得到上三角系统 $U\vec{x}=\vec{y}$ 后，从最后一行开始求解：

BackSubstitution(U, y):
    # x 先保存右端项，随后逐步改写为解。
    x <- y

    for p <- n, n - 1, ..., 1:
        # 先解出当前主元对应的未知量。
        x_p <- x_p / u_pp

        # 再消去第 p 列在主元上方的元素。
        for r <- 1, 2, ..., p - 1:
            x_r <- x_r - u_rp * x_p

    return x

复杂度：$\text{Gaussian elimination}=O(n^3)$，$\text{triangular solve}=O(n^2)$。

这也是 LU 分解有用的原因：同一个系数矩阵 $A$ 若需要搭配多个右端项 $\vec{b}$ 反复求解，可以先花一次 $O(n^3)$ 做分解，之后每次只需 $O(n^2)$。

LU 分解

设一系列消元矩阵把 $A$ 化为上三角矩阵：

$$M_k\cdots M_1A=U$$

则：

$$A=(M_k\cdots M_1)^{-1}U$$

由：

$$(AB)^{-1}=B^{-1}A^{-1}$$

得到：

$$A=(M_1^{-1}M_2^{-1}\cdots M_k^{-1})U$$

定义：

$$L=M_1^{-1}M_2^{-1}\cdots M_k^{-1}$$

于是：

$$A=LU$$

消元矩阵的逆是下三角矩阵，下三角矩阵乘积仍是下三角矩阵，因此 $L$ 是下三角矩阵；消元后的 $U$ 是上三角矩阵。

使用 LU 求解

若 $A=LU$，则：

$$A\vec{x}=\vec{b}\quad\Longleftrightarrow\quad LU\vec{x}=\vec{b}$$

令 $\vec{y}=U\vec{x}$。

求解分两步：

1. 先解下三角系统：

$$L\vec{y}=\vec{b}$$

2. 再解上三角系统：

$$U\vec{x}=\vec{y}$$

实际计算中通常不显式构造 $L^{-1}$ 或 $U^{-1}$，而是用前向代入和后向代入。

置换与 LU

不是所有矩阵都能在不置换的情况下稳定写成 $A=LU$。若消元需要换行，更常见的形式是：

$$PA=LU$$

或等价地：

$$A=P^TLU$$

其中 $P$ 是置换矩阵。

紧凑存储

若约定 $L$ 的对角线元素全部为 1，则不必显式存储 $L$ 的对角线。$L$ 的严格下三角部分和 $U$ 的上三角部分可以共享原矩阵空间：

$$\begin{bmatrix} U & U & U & U\\ L & U & U & U\\ L & L & U & U\\ L & L & L & U\\ L & L & L & L \end{bmatrix}$$

该图表示存储位置归属：上三角区域存放 $U$，下三角区域存放 $L$ 的非对角元素，$L$ 的对角线默认为 1。

线性系统建模

线性系统不仅用于直接求解方程，也可以用于从数据中拟合模型。回归问题的目标是根据实验输入 $\vec{x}$ 和观测输出 $y$，构造函数 $f$，用于预测新输入的输出。

线性回归

最简单的线性模型为：

$$f(\vec{x})=a_1x_1+a_2x_2+\cdots+a_nx_n=\vec{a}^{\,T}\vec{x}$$

若有多组实验数据 $\{(\vec{x}^{(i)},y^{(i)})\}_{i=1}^{m}$，则可以写成线性系统：

$$A\vec{a}=\vec{y}$$

其中 $A$ 的第 $i$ 行是 $\vec{x}^{(i)T}$。当 $m=n$ 且 $A$ 可逆时，可以精确解出参数；当 $m>n$ 时，一般转化为最小二乘问题。

更一般地，模型可以用基函数线性组合表示：

$$f(\vec{x})=a_1f_1(\vec{x})+a_2f_2(\vec{x})+\cdots+a_nf_n(\vec{x})$$

虽然 $f_i$ 本身可以是非线性的，但未知参数 $a_i$ 仍然线性出现，因此仍可写成线性系统。

多项式回归

一维多项式回归使用基函数 $1,x,x^2,\dots,x^{n-1}$：

$$f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}$$

对应的系数矩阵是 Vandermonde 矩阵：

$$\begin{bmatrix} 1 & x^{(1)} & (x^{(1)})^2 & \cdots & (x^{(1)})^{n-1}\\ 1 & x^{(2)} & (x^{(2)})^2 & \cdots & (x^{(2)})^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x^{(m)} & (x^{(m)})^2 & \cdots & (x^{(m)})^{n-1} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix} = \begin{bmatrix} y^{(1)}\\ y^{(2)}\\ \vdots\\ y^{(m)} \end{bmatrix}$$

精确插值会强制模型通过所有数据点。如果观测含噪声，或者基函数选得不合适，精确拟合可能导致过拟合。

最小二乘与正规方程

当 $m>n$ 时，$A\vec{x}=\vec{b}$ 通常不可精确满足。最小二乘选择残差平方最小的解：

$$\min_{\vec{x}}\|A\vec{x}-\vec{b}\|_2^2$$

展开目标函数：

$$\|A\vec{x}-\vec{b}\|_2^2=\vec{x}^{\,T}A^TA\vec{x}-2(A^T\vec{b})\cdot\vec{x}+\|\vec{b}\|_2^2$$

令梯度为 0，得到正规方程：

$$A^TA\vec{x}=A^T\vec{b}$$

这个方程将一个通常无精确解的超定线形方程转化成一个可以求解的方程组。其中 $A^TA$ 称为 Gram 矩阵。若 $A$ 的列线性无关，则 $A^TA$ 正定，正规方程有唯一解。

正则化

欠定问题或病态问题需要额外约束。Tikhonov 正则化在残差项外加入参数范数惩罚：

$$\min_{\vec{x}}\|A\vec{x}-\vec{b}\|_2^2+\alpha\|\vec{x}\|_2^2,\qquad 0<\alpha\ll 1$$

其正规方程为：

$$(A^TA+\alpha I)\vec{x}=A^T\vec{b}$$

$\alpha$ 增大会改善可逆性和数值稳定性，但解不再精确满足 $A\vec{x}=\vec{b}$，并且会引入偏差。

常见正则化形式：

• Tikhonov：$\min_{\vec{x}}\|A\vec{x}-\vec{b}\|_2^2+\alpha\|\vec{x}\|_2^2$。
• Lasso：$\min_{\vec{x}}\|A\vec{x}-\vec{b}\|_2^2+\beta\|\vec{x}\|_1$。
• Elastic Net：$\min_{\vec{x}}\|A\vec{x}-\vec{b}\|_2^2+\alpha\|\vec{x}\|_2^2+\beta\|\vec{x}\|_1$。

Gram 矩阵与 Cholesky 分解

Gram 矩阵 $A^TA$ 总是对称半正定，因为：

$$\vec{x}^{\,T}(A^TA)\vec{x}=(A\vec{x})^T(A\vec{x})=\|A\vec{x}\|_2^2\ge 0$$

而当 $A^TA$ 正定时，由 $\|A\vec{x}\|_2^2 > 0$ 有 $\forall x\ne 0, Ax \ne 0$。这个结论则等价于“矩阵 $A$ 的所有列向量线性无关”，即 $A$ 列满秩。

基于以上得出的结论，可以引入 Cholesky 分解方法。若 $C$ 对称正定，则可以分解为：

$$C=LL^T$$

其中 $L$ 是下三角矩阵。Cholesky 分解相当于专门针对对称正定矩阵的 LU 分解，只需存储一半矩阵，它就像正定矩阵的一种开平方。

假设下三角矩阵 $L$ 为：

$$L = \begin{bmatrix}l_{11} \\ l_{21} & l_{22} \\ l_{31} & l_{32} & l_{33}\end{bmatrix}$$

那么有：

$$LL^T = \begin{bmatrix}l_{11}^2 & l_{11}l_{21} & l_{11}l_{31} \\ l_{21}l_{11} & l_{21}^2 + l_{22}^2 & l_{21}l_{31} + l_{22}l_{32} \\ l_{31}l_{11} & l_{31}l_{21} + l_{32}l_{22} & l_{31}^2 + l_{32}^2 + l_{33}^2\end{bmatrix}$$

按顺序从左到右从上到下列等式计算即可，写成通式形式如下：

1. 先计算第 $k$ 个对角元：

   $$L_{kk} = \sqrt{C_{kk} - \sum_{j=1}^{k-1}L_{kj}^2}$$

2. 然后计算第 $k$ 列下面的元素：

   $$L_ik = \frac{C_{ik}-\sum_{j=1}^{k-1}L_{ij}L_{kj}}{L_{kk}} \quad i = k+1, ..., n$$

写代码这么写，手算不如直接推一遍。

算法流程：

Cholesky(C):
    # C 对称正定，目标是 C = L L^T。
    L <- 0

    for k <- 1, 2, ..., n:
        # 对角元由当前剩余平方根给出。
        L_kk <- sqrt(C_kk - sum_{j<k} L_kj^2)

        # 第 k 列下方元素由对应内积修正后除以 L_kk。
        for i <- k + 1, ..., n:
            L_ik <- (C_ik - sum_{j<k} L_ij L_kj) / L_kk

    return L

范数与条件数

线性系统中的“误差很小”必须依赖范数定义。范数用于度量向量或矩阵大小。

向量范数

向量范数 $\|\cdot\|:\mathbb{R}^n\to[0,\infty)$ 满足：

• 非负性：$\|\vec{x}\|=0$ 当且仅当 $\vec{x}=\vec{0}$。
• 齐次性：$\|c\vec{x}\|=|c|\|\vec{x}\|$。
• 三角不等式：$\|\vec{x}+\vec{y}\|\le \|\vec{x}\|+\|\vec{y}\|$。

常见 $p$-范数：

$$\|\vec{x}\|_p=\left(|x_1|^p+|x_2|^p+\cdots+|x_n|^p\right)^{1/p},\qquad p\ge 1$$

特殊情形包括：

$$\|\vec{x}\|_2=\sqrt{\sum_i x_i^2},\qquad \|\vec{x}\|_1=\sum_i|x_i|,\qquad \|\vec{x}\|_\infty=\max_i|x_i|$$

有限维空间中所有范数等价，即不同范数给出的“小”和“大”只差常数因子，选取不同的范数不会改变对矩阵“大”和“小”的判断。

矩阵范数

一种构造方式是把矩阵展开。例如 Frobenius 范数，其直接计算所有权重绝对值平方总和的开方：

$$\|A\|_{\mathrm{Fro}}=\sqrt{\sum_{i,j}a_{ij}^2}$$

另一种更重要的构造是诱导范数。矩阵在几何层面上能对图形进行拉伸，这种构造方式利用矩阵最大能将一个单位向量拉到多长来判断其“大小”：

$$\|A\|=\max_{\|\vec{x}\|=1}\|A\vec{x}\|$$

一范数对每一列绝对值求和，然后取最大值；无穷范数对每一行绝对值求和，然后取最大值：

$$\|A\|_1=\max_j\sum_i|a_{ij}|,\qquad \|A\|_\infty=\max_i\sum_j|a_{ij}|$$

对 2-范数，$\|A\|_2$ 等于最大奇异值，也就是单位球被 $A$ 映射后的最大拉伸长度。

矩阵条件数

对可逆矩阵 $A$，矩阵条件数定义为：

$$\operatorname{cond}(A)=\kappa(A)=\|A\|\,\|A^{-1}\|$$

若 $A$ 不可逆，则 $\operatorname{cond}(A)=\infty$，因为某些非 0 方向会被压成 0，线性系统可能没有唯一解，对扰动及其敏感。条件数衡量线性系统对输入扰动的敏感性。

考虑扰动问题：

$$(A+\varepsilon\delta A)\vec{x}(\varepsilon)=\vec{b}+\varepsilon\delta\vec{b}$$

令相对扰动量：

$$D=\frac{\|\delta\vec{b}\|}{\|\vec{b}\|}+\frac{\|\delta A\|}{\|A\|}$$

则解的相对变化满足一阶估计：

$$\frac{\|\vec{x}(\varepsilon)-\vec{x}(0)\|}{\|\vec{x}(0)\|}\le \varepsilon D\kappa(A)+O(\varepsilon^2)$$

条件数越大，同样的输入扰动越可能导致更大的解误差。对诱导范数，条件数也可以理解为最大拉伸与最小拉伸之比：

$$\operatorname{cond}(A)=\frac{\max_{\|\vec{x}\|=1}\|A\vec{x}\|}{\min_{\|\vec{x}\|=1}\|A\vec{x}\|}$$

这说明条件数衡量的是矩阵在不同方向上拉伸程度的不均匀性。

列空间与 QR 分解

最小二乘问题的几何解释是：把 $\vec{b}$ 投影到 $A$ 的列空间中，找到最接近 $\vec{b}$ 的向量 $A\vec{x}$。

标准正交矩阵

若 $Q$ 的列向量标准正交，则：

$$Q^TQ=I$$

因此 $Q$ 保持内积与长度：

$$(Q\vec{x})\cdot(Q\vec{y})=\vec{x}\cdot\vec{y},\qquad \|Q\vec{x}\|_2=\|\vec{x}\|_2$$

正交矩阵适合数值计算，因为它不会放大向量长度，也不会恶化条件数。

QR 分解

QR 分解的目的在于把矩阵 $A$ 的列空间换成一组更好用的标准正交基 $Q$，再用上三角矩阵 $R$ 记录原来每一列在这些正交基下的坐标。

QR 分解把矩阵写成：

$$A=QR$$

其中 $Q$ 的列标准正交，$R$ 是上三角矩阵。若 $A=QR$，最小二乘正规方程变为：

$$A^TA\vec{x}=A^T\vec{b}\quad\Longleftrightarrow\quad R^TR\vec{x}=R^TQ^T\vec{b}$$

当 $R$ 可逆时，等价于：

$$\vec{x}=R^{-1}Q^T\vec{b}$$

如果直接解正规方程，需要显式构造 $A^TA$ ，其条件数为 $A$ 条件数的平方，如果 $A$ 本来就病态特征，会被 $A^TA$ 进一步放大。这样避免显式构造 $A^TA$，通常比正规方程更稳定。

Gram-Schmidt 正交化

向量 $\vec{b}$ 在非零向量 $\vec{a}$ 方向上的投影为：

$$\operatorname{proj}_{\vec{a}}\vec{b}=\frac{\vec{a}^{\,T}\vec{b}}{\|\vec{a}\|_2^2}\vec{a}$$

残差 $\vec{b}-\operatorname{proj}_{\vec{a}}\vec{b}$ 与 $\vec{a}$ 正交。这是最小二乘投影的基本几何事实。

若 $\hat{a}_1,\dots,\hat{a}_k$ 标准正交，则投影到其张成空间为：

$$\operatorname{proj}_{\operatorname{span}\{\hat{a}_1,\dots,\hat{a}_k\}}\vec{b}=\sum_{i=1}^{k}(\hat{a}_i^{\,T}\vec{b})\hat{a}_i$$

给定向量 $\vec{v}_1,\dots,\vec{v}_k$，Gram-Schmidt 依次去掉当前向量在已有正交基方向上的投影：

$$\hat{a}_1=\frac{\vec{v}_1}{\|\vec{v}_1\|}$$

对 $i=2,\dots,k$：

$$\vec{p}_i=\operatorname{proj}_{\operatorname{span}\{\hat{a}_1,\dots,\hat{a}_{i-1}\}}\vec{v}_i,\qquad \hat{a}_i=\frac{\vec{v}_i-\vec{p}_i}{\|\vec{v}_i-\vec{p}_i\|}$$

按照施密特正交化的顺序，原 $A$ 矩阵的第 $k$ 列只由 $Q$ 矩阵的前 $k$ 列线性组合而成。可以推出在矩阵 $R$ 中，有：

$$r_{ij} = \hat{a}_i^Ta_j, \quad i < j$$

$$r_{ii} = ||u_i|| = ||v_i - \sum_{j=1}^{i-1}(\hat{a}_j^Tv_i)\hat{a}_j||$$

算法流程：

GramSchmidt(v_1, ..., v_k):
    for i <- 1, 2, ..., k:
        u_i <- v_i
        
        # 去掉在已有正交方向上的分量。
        for j <- 1, 2, ..., i - 1:
            u_i <- u_i - (a_j^T v_i) a_j
            
        # 归一化得到新的正交单位向量。
        a_i <- u_i / ||u_i||
        
    return a_1, ..., a_k

经典 Gram-Schmidt 容易受数值误差影响。改进 Gram-Schmidt 会在每一步用更新后的残差继续投影，通常更稳定。

Householder QR

Householder 反射是一类容易构造的正交矩阵：

$$H_{\vec{v}}=I-\frac{2\vec{v}\vec{v}^{\,T}}{\vec{v}^{\,T}\vec{v}}$$

它可以把向量反射到指定坐标轴方向，从而一次消去某一列对角线下方的所有元素。该算法的关键在于逐列三角化，Householder QR 的基本流程：

HouseholderQR(A):
    Q <- I
    R <- A

    for k <- 1, 2, ..., n:
        # 取出 R 的第 k 列从第 k 行开始的部分。
        a <- tail of column k

        # 构造反射向量 v，使 H_v a 只保留第一个分量。
        v <- a - ||a|| e_1

        # 左乘反射矩阵，消去第 k 列对角线下方元素。
        R <- H_v R
        Q <- Q H_v^T

    return Q, R

实际实现中通常不显式存储完整 $Q$，而是存储反射向量 $\vec{v}$。

简化 QR

若 $A\in\mathbb{R}^{m\times n}$ 且 $m\gg n$，完整 QR 中 $Q\in\mathbb{R}^{m\times m}$ 过大。常用 reduced QR：

$$A=Q_1R_1,\qquad Q_1\in\mathbb{R}^{m\times n},\quad R_1\in\mathbb{R}^{n\times n}$$

只保留与列空间相关的正交基即可。

特征值与特征向量

特征值问题研究矩阵在某些方向上的纯缩放作用。

基本定义

若 $\vec{x}\ne\vec{0}$ 且：

$$A\vec{x}=\lambda\vec{x}$$

则 $\lambda$ 是特征值，$\vec{x}$ 是对应特征向量。复数域中，每个 $A\in\mathbb{C}^{n\times n}$ 至少有一个特征值。非 0 向量 $\vec{x}$ 在 $A$ 的变换下不会改变方向，只可能出现长度的拉伸，$\lambda$ 指示了向量长度放缩的倍数，若特征值为负则说明特征向量经过变换反向。特征值问题研究矩阵在哪些方向上表现为一个纯缩放因子。

按照一下步骤解特征值：

$$Ax = \lambda x$$

$$(A - \lambda I)x = 0 \quad x \ne 0$$

齐次方程有唯一解当且仅当：

$$det(A - \lambda I) = 0$$

根据该特征方程解特征值。

谱 $\sigma(A)$ 是全部特征值的集合；谱半径为：

$$\rho(A)=\max_i|\lambda_i|$$

若矩阵有 $n$ 个线性无关特征向量，则它可对角化：

$$A=XDX^{-1}$$

其中 $D$ 是特征值对角矩阵，$X$ 的列为特征向量。于是：

$$A^{-1}=XD^{-1}X^{-1}$$

将 $A$ 分解得到 $D$ 的过程即对角化。在特征向量组成的坐标系中，矩阵 $A$ 的作用只是每个方向独立缩放。对角化后幂和逆都很容易计算。

Hermitian 矩阵与谱定理

复矩阵的共轭转置记为 $A^H$。若：

$$A=A^H$$

则 $A$ 是 Hermitian 矩阵。实数情形下就是对称矩阵。

$$x^H A x = \lambda x^H x$$

$x^H A x$ 一定等于自己的复共轭，其一定是实数，又 $x^H x = ||x||^2 > 0$，$\lambda = \frac{x^H A x}{x^H x}$ 一定是实数。

Hermitian 矩阵具有重要性质：

• 所有特征值都是实数。
• 不同特征值对应的特征向量相互正交。
• 存在一组标准正交特征向量。

谱定理：

• 实数形式下，任意实对称矩阵有分解：

  $$A=XDX^T$$

  其中 $A$ 是实对称矩阵， $X$ 是正交矩阵，$D$ 是实对角矩阵。

• 复数形式下，对于满足 $A = A^H$ 的厄米矩阵 $A$ 有分解：

  $$A = U D U^H$$

  其中 $U$ 的列向量是一组标准正交特征向量，$D$ 的对角线元素是实特征值。

Rayleigh 商

对于对称矩阵，最大化 $\vec{x}^{\,T}A\vec{x}$ 且约束 $\|\vec{x}\|_2=1$，可得到特征值问题：

$$A\vec{x}=\lambda\vec{x}$$

对应目标值为：

$$\vec{x}^{\,T}A\vec{x}=\lambda$$

Rayleigh 商定义为：

$$R_A(\vec{x})=\frac{\vec{x}^{\,T}A\vec{x}}{\vec{x}^{\,T}\vec{x}}$$

它常用于估计给定近似特征向量对应的特征值。可以理解为：给定一个向量 $\vec{x}$，估计它“像哪个特征值对应的特征向量”。

幂迭代

设 $A$ 对称，特征值满足 $|\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_n|$，初始向量 $\vec{v}_0$ 在主特征向量方向上的分量非零。幂迭代：

$$\vec{w}_k=A\vec{v}_{k-1},\qquad \vec{v}_k=\frac{\vec{w}_k}{\|\vec{w}_k\|}$$

会收敛到最大模特征值对应的特征向量方向。其直观原因是：

$$A^k\vec{v}\approx \lambda_1^k c_1\vec{x}_1$$

要证明其收敛性，展开左式，将 $\vec{v}$ 张开成基向量的线性组合：

$$A^k v_0 = c_1\lambda_1^k x_1 + c_2\lambda_2^k x_2 + \dots + c_n\lambda_n^k x_n$$

提出 $\lambda_1$：

$$A^k v_0 = \lambda_1^k \left[ c_1x_1 + c_2\left(\frac{\lambda_2}{\lambda_1}\right)^k x_2 + \dots + c_n\left(\frac{\lambda_n}{\lambda_1}\right)^k x_n \right]$$

随着 $k$ 的增大，第一项意外的所有项会向 0 收敛，最终稳定到 $\lambda_1^k c_1\vec{x}_1$。由于除去 $\lambda_1$ 之外最大模特征值为 $\lambda_2$，这个收敛速度就由 $|\lambda_2/\lambda_1|$ 来控制，当 $\lambda_2$ 对应的项收敛到 0 附近时，后面的项也一定收敛至 0 附近了。

反迭代与位移反迭代

若 $A\vec{v}=\lambda\vec{v}$，则：

$$A^{-1}\vec{v}=\frac{1}{\lambda}\vec{v}$$

因此对 $A^{-1}$ 做幂迭代可以找到 $A$ 的最小模特征值对应方向。实际不显式求逆，而是每步解线性系统：

InverseIteration(A, v_0):
    factor A as LU

    for k <- 1, 2, ...:
        solve A w_k = v_{k-1}
        v_k <- w_k / ||w_k||

    return v_k

若想找到最接近 $\sigma$ 的特征值，可使用位移反迭代：

$$\vec{w}_k=(A-\sigma I)^{-1}\vec{v}_{k-1},\qquad \vec{v}_k=\frac{\vec{w}_k}{\|\vec{w}_k\|}$$

因为 $A-\sigma I$ 的特征值为 $\lambda_i-\sigma$，离 $\sigma$ 最近的 $\lambda_i$ 会在逆矩阵中变成最大模特征值。

Rayleigh 商迭代每步更新位移：

$$\sigma_{k+1}=\frac{\vec{v}_k^{\,T}A\vec{v}_k}{\|\vec{v}_k\|_2^2}$$

并求解：

$$\vec{w}_k=(A-\sigma_kI)^{-1}\vec{v}_{k-1},\qquad \vec{v}_k=\frac{\vec{w}_k}{\|\vec{w}_k\|}$$

奇异值分解

特征值分解解决了方阵的特征分解问题，但是这种方法不适用于非方阵的矩阵。SVD 描述一般矩阵的几何作用：先在输入空间旋转，再沿坐标轴缩放，最后在输出空间旋转。换言之，我们要在输入空间找一组标准正交方向，在输出空间找一组标准正交方向，使得 $A$ 只是在这些方向之间做缩放。

从特征值到奇异值

对任意 $A\in\mathbb{R}^{m\times n}$，矩阵 $A^TA$ 对称半正定。若：

$$A^T A\vec{v}_i=\lambda_i\vec{v}_i,\qquad \lambda_i\ge 0$$

则奇异值定义为：

$$\sigma_i=\sqrt{\lambda_i}$$

即 $A^TA$ 的特征值。当 $\sigma_i>0$ 时，对应左奇异向量为：

$$\vec{u}_i=\frac{A\vec{v}_i}{\sigma_i}$$

改写成 $AV = U\Sigma$ 的形式表示：矩阵 $A$ 把输入空间中的方向 $v_i$​，映射到输出空间中的方向 $u_i$​，并且长度放大 $\sigma_i$ 倍。将右奇异向量组成 $V$，左奇异向量组成 $U$，奇异值组成对角矩阵 $\Sigma$，得到：

$$A=U\Sigma V^T$$

其中 $U$ 和 $V$ 正交，$\Sigma$ 的对角线为 $\sigma_1\ge\sigma_2\ge\cdots\ge 0$。

SVD 的常用解释：

• $V$ 的列是右奇异向量，张成输入空间的正交基。
• $U$ 的列是左奇异向量，张成输出空间的正交基。
• $\Sigma$ 表示沿各主方向的缩放。
• 非零奇异值对应 $A$ 的列空间与行空间。

用 SVD 解线性系统

若 $A=U\Sigma V^T$，则：

$$A\vec{x}=\vec{b}\quad\Longleftrightarrow\quad \Sigma\vec{y}=U^T\vec{b},\qquad \vec{x}=V\vec{y}$$

这样通过 SVD 将原来的超定方程变换成了一个对角系统。$d = \Sigma y$ 对每个奇异值：

$$y_i=\begin{cases} d_i/\sigma_i, & \sigma_i\ne 0\\ 0, & \sigma_i=0 \end{cases}$$

其中 $\vec{d}=U^T\vec{b}$。这自然给出伪逆。

伪逆

若：

$$A=U\Sigma V^T$$

则 Moore-Penrose 伪逆为：

$$A^+=V\Sigma^+U^T$$

其中：

$$\Sigma^+_{ij}=\begin{cases} 1/\sigma_i, & i=j,\ \sigma_i\ne 0\\ 0, & \text{otherwise} \end{cases}$$

伪逆的作用：

• 若 $A$ 方阵可逆，则 $A^+=A^{-1}$。
• 若 $A$ 高矩阵，$A^+\vec{b}$ 给出最小二乘解。
• 若 $A$ 宽矩阵，$A^+\vec{b}$ 给出所有最小二乘解中欧几里得范数最小的解。

SVD 与矩阵范数

SVD 给出常用矩阵范数的简单表达：

$$\|A\|_{\mathrm{Fro}}^2=\sum_i\sigma_i^2$$

$$\|A\|_2=\max_i\sigma_i$$

若 $A$ 满秩且可逆，则 2-范数条件数为：

$$\operatorname{cond}_2(A)=\frac{\sigma_{\max}}{\sigma_{\min}}$$

很小的奇异值会导致伪逆中出现很大的 $1/\sigma_i$，从而放大噪声。实际计算中常忽略过小奇异值来稳定结果。

低秩近似

SVD 可以写成外积和：

$$A=\sum_{i=1}^{r}\sigma_i\vec{u}_i\vec{v}_i^{\,T}$$

截断到前 $k$ 个最大奇异值：

$$A_k=\sum_{i=1}^{k}\sigma_i\vec{u}_i\vec{v}_i^{\,T}$$

Eckart-Young 定理说明，$A_k$ 是秩不超过 $k$ 的最佳近似，同时最小化 Frobenius 范数误差和 2-范数误差：

$$A_k=\arg\min_{\operatorname{rank}(B)\le k}\|A-B\|_{\mathrm{Fro}}$$

SVD 与正则化

Tikhonov 正则化问题：

$$(A^TA+\alpha I)\vec{x}=A^T\vec{b}$$

若 $A=U\Sigma V^T$，则解可以写成：

$$\vec{x}=V D U^T\vec{b},\qquad D_{ii}=\frac{\sigma_i}{\sigma_i^2+\alpha}$$

正则化会把原本的 $1/\sigma_i$ 替换为 $\sigma_i/(\sigma_i^2+\alpha)$，从而抑制小奇异值带来的噪声放大。

PCA

PCA 可以表述为寻找正交投影子空间，使数据重构误差最小。若数据矩阵为 $X$，目标是：

$$\min_{C^TC=I}\|X-CC^TX\|_{\mathrm{Fro}}$$

若 $X=U\Sigma V^T$，最优的 $C$ 由 $U$ 的前 $d$ 列给出。也就是说，PCA 的主方向对应最大奇异值方向。