矩阵范数

Gauss 消元

    Gauss 消元 设 `bm A` 是 `n` 阶可逆实方阵, `bm b` 是 `n` 维列向量. 用如下步骤求解线性方程组 `bm (A x) = bm b`:
  1. 将矩阵 `bm A` 与 `bm b` 合成一个增广矩阵 `[bm A, bm b]`.
  2. 从 `bm A` 的第一列找到绝对值最大的元素 `p = a_(i 1)`, 称为列主元. 由于 `bm A` 可逆, `|p| gt 0`.
  3. 将增广矩阵的第一行与第 `i` 行交换, 使得列主元占据左上角的位置.
  4. 对第一列进行 Gauss 消元. 将增广矩阵第一行的 `-a_(i 1)//p` 倍加到第 `i` 行, `i = 2, cdots, n`, 使得矩阵第一列除第一行外全为零, 形如 `[p, ast, ast; 0, bm A_2, bm b_2]`.
  5. 继续对右下角的矩阵 `[bm A_2, bm b_2]` 使用前面的步骤, 选主元并消元, 最终矩阵化为 `[bm R, bm c]`, 其中 `bm R` 是上三角矩阵, 对角元非零.
  6. 使用代入法轻松求解上三角形的方程组 `bm (R x) = bm c`: `x_n = c_n//r_(n n)`,
    `x_i = (c_i - sum_(j=i+1)^n r_(i j) x_j) // r_(i i)`, `quad i = n-1, cdots, 1`.

求矩阵的秩 模仿 Gauss 消元的流程, 但只需化为上三角矩阵即可, 不需要代入的步骤.

    一般消元法 对任意实矩阵 `bm A_(m xx n)` 和任意 `n` 维向量 `bm b` 求解线性方程组 `bm (A x) = bm b`.
  1. 选主元. 和 Gauss 消元一样, 从增广矩阵 `[bm A, bm b]` 开始, 遍历 `bm A` 的各列, 选出当前列的绝对值最大者作为列主元 (跳过那些已被选为主元的行).
  2. 消元步骤. 记第 `j` 列的主元是 `p = a_(i j)`. 首先将第 `i` 行同乘以 `1//p` 使主元化为 1, 再将这一行的 `-a_(k j)` 倍加到第 `k` 行, `k = 1, cdots, i-1, i+1, cdots, m`. 这样操作后, 第 `j` 列除主元所在行等于 1 以外, 其余元素均为 0.
  3. 退化情形. 当然, 有可能某一列的元素全为 0, 或者矩阵的每一行都已经选出主元, 导致当前列不存在有效主元. 在这两种情况下, 不能进行 2. 的消元步骤. 如果某一列选不出有效主元, 就称这一列对应的未知数为非基变量, 否则称为基变量.
  4. 解的存在. 理想情况下, 若 `bm A` 行满秩, 则每一行都能选出一个主元. 否则, 需要检查未能选出主元的那些行, 如果该行的常数项不为 0, 则方程无解.
  5. 求解步骤. 若方程有解, 把基变量留在方程左边, 非基变量移动到方程右边. 每个基变量的系数恰好都是 1, 方程形如 `bm (I x_1) = bm c - bm (B x_2)`, 其中 `bm x_1, bm x_2` 分别为基变量和非基变量, `bm c` 是常数项, `bm B` 是非基变量的系数矩阵. 方程的解已经显化出来: `[bm x_1; bm x_2] = [bm c - bm B bm x_2; bm x_2]`, `quad` 其中 `bm x_2` 取遍 `RR^(n-r)`.

注: 输出的第一列是特解, 后面各列是基础解系.

QR 分解

QR 分解 [来自 知乎@Iterator] 设 `bm A` 是 `m xx n` 实矩阵, 且 `bm A` 的各列向量线性无关, 则存在唯一分解 `bm A_(m xx n) = bm Q_(m xx n) bm R_(n xx n)`. 其中 `bm Q` 满足各列向量单位正交, `bm R` 是主对角线全为正数的上三角矩阵. QR 分解以矩阵的形式, 编码了 Schmidt 正交化的计算过程.

  1. 存在性. 对 `bm A` 的各列进行 Schmidt 正交化, 其中 `k_(i,j) = ((bm a_i, bm b_j))/((bm b_j, bm b_j))`: `bm b_1 = bm a_1`,
    `bm b_2 = bm a_2 - k_(2,1) bm b_1`,
    `cdots`
    `bm b_n = bm a_n - k_(n,1) bm b_1 - k_(n,2) bm b_2 - cdots - k_(n,n-1) bm b_(n-1)`.     再单位化, `hat bm b_j = bm b_j/|bm b_j|`, `quad j = 1, cdots, n`. 写为矩阵形式就是 `bm A = (hat bm b_1, cdots, hat bm b_n)[ |bm b_1|, k_(2,1)|bm b_1|, cdots, k_(n,1)|bm b_1|; , |bm b_2|, cdots, k_(n,2)|bm b_2|; ,,ddots,vdots; ,,,|bm b_n| ]`.
  2. 唯一性. 设 `bm A` 有 QR 分解 `bm A = bm Q_1 bm R_1 = bm Q_2 bm R_2`, 下证 `bm R_1 = bm R_2`.
    由于 `bm Q_1, bm Q_2` 的各列向量单位正交, 有 `bm Q_1^T bm Q_1 = bm Q_2^T bm Q_2 = bm I`, 于是 `bm R_1^T bm R_1` `= bm R_1^T bm Q_1^T bm Q_1 bm R_1` `= bm A^T bm A` `= bm R_2^T bm Q_2^T bm Q_2 bm R_2` `= bm R_2^T bm R_2`. 上式两边左乘 `(bm R_2^T)^-1`, 右乘 `bm R_1^-1` 得 `(bm R_2^T)^-1 bm R_1^T = bm R_2 bm R_1^-1`. 上式左边是下三角矩阵, 右边是上三角矩阵, 但两边相等, 从而只能是对角形矩阵, 记为 `bm D`. 则 `bm R_2 = bm D bm R_1`. 于是 `bm R_1^T bm R_1` `= bm R_2^T bm R_2` `= bm R_1^T bm D^T bm D bm R_1`, 即 `bm D^T bm D = bm I`. 显然 `bm D` 的主对角线全为正, 故 `bm D = bm I`, 即 `bm R_1 = bm R_2`. 从而 `bm Q_1 = bm A bm R_1^-1 = bm A bm R_2^-1 = bm Q_2`.

QR 分解的 Q 源于 orthogonal 的 O, 为避免混淆而写作 Q; 而 R 源于 right triangular matrices. 

# sympy 求 QR 分解
A.QRdecomposition()

注: 若 `bm A` 的列向量线性相关, 则 `bm Q` 含有零列, `bm R` 含有零行.

在最小二乘法中, 对于可能无解的实系数方程组 `bm (A x) = bm b`, 两边左乘 `bm A^T`, 转化为 `bm A^T bm A x = bm A^T bm b`, 我们需要计算 `bm x = (bm A^T bm A)^-1 bm A^T bm b`. 如果先将 `bm A` 分解为 `bm A = bm (Q R)`, 则上式化为 `bm x` `= (bm R^T bm Q^T bm Q bm R)^-1 bm R^T bm Q^T bm b` `= (bm R^T bm R)^-1 bm R^T bm Q^T bm b` `= bm R^-1 bm Q^T bm b`. 从而将复杂度较高的矩阵求逆运算, 化为较容易的上三角矩阵 `bm R` 的求逆运算.

QR 分解 (基于 Householder 变换) Householder 变换是指镜面反射, 或称反向正交变换. 它可以写为 `bm H(bm u) = bm I - 2 bm u bm u^T` 的形式, 其中 `bm u` 为单位向量. 由这里知道, 已知向量 `bm x` 和单位向量 `bm y`, 要使 `bm (H x)` 平行于 `bm y`, 只需取 `bm H = bm H(hat bm u)`, 其中 `bm u = bm x - |bm x| bm y`, `quad hat bm u = bm u //|bm u|`. 记实矩阵 `bm A` 的第 `j` 列为 `bm a_j`, 又记 `bm e_i` 是第 `i` 个分量为 1, 其余分量为 0 的单位向量. 先作 Householder 变换 `bm H_1` 使得 `bm (H_1 bm a_1)` 平行于 `bm e_1`, 得到 `bm (H_1 A)` `= [|bm a_1|, ast; 0, bm A_2]`. 对于右下角的 `bm A_2`, 继续以 Householder 变换 `bm H_2` 将它的第一列除第一行外的所有元素化为零. 最终得到上三角矩阵 `bm R = bm H_n cdots bm H_1 bm A`, 其中 `bm Q = bm H_n cdots bm H_1` 是正交矩阵.

特征值问题

下面这个定理可以帮助我们估计特征值的范围:

Gerschgorin (盖尔) 圆盘定理 设 `bm A in CC^(n xx n)`, 分别定义 行盖尔圆盘: `S_i := {z in CC: |z-a_(i i)| le sum_(j!=i) |a_(i j)|}`,
列盖尔圆盘: `G_i := {z in CC: |z-a_(i i)| le sum_(j!=i) |a_(j i)|}`. 行 (列) 盖尔圆盘的半径分别叫做去心行 (列) 和. `bm A` 的任一特征值 `lambda` 必然落在某个行盖尔圆盘中, 也必然落在某个列盖尔圆盘中: `lambda in (uuu_(i=1)^n S_i) nn (uuu_(i=1)^n G_i)`.

以行盖尔圆盘为例, 任取 `bm A` 的特征值 `lambda` 和特征向量 `bm x != 0`, 又设 `x_i` 是 `bm x` 的各分量中模最大的一个, 则由 `bm (A x) = lambda bm x` 知, `sum_(j=1)^n a_(i j) x_j = lambda x_i`, 于是 `|x_i(a_(i i) - lambda)|` `= |sum_(j!=i) a_(i j) x_j|` `le |x_i| sum_(j!=i) |a_(i j)|`. 这说明 `lambda` 落在一个行盖尔圆盘中.

    Jacobi 迭代法 适用于求实对称矩阵 `bm A` 的特征值与特征向量. 思想是作相似变换 `bm T^-1 bm A bm T`, 使得变换后的矩阵接近对角矩阵. 其中, `bm T = bm T(i, j, theta)` 是一个 Givens 旋转矩阵, 满足 `[T_(i i), T_(i j); T_(j i), T_(j j)]` `= [cos theta, -sin theta; sin theta, cos theta]`, 其它元素与单位矩阵 `bm I` 相同. `i, j` 元的选择应当使得 `|A_(i j)|` 在非对角元上取得最大值, 称为主元 (pivot). `theta` 的选择应当使得 `bm A^((1)) = bm T^-1 bm A bm T` 的 `i j` 元等于零, 即 `A^((1))_(i j) = A_(i j) cos 2 theta + (A_(i i) - A_(j j))/2 sin 2 theta = 0`. 于是 `theta = {1/2 arctan((2 A_(i j))/(A_(i i) - A_(j j))), if A_(i i) != A_(j j); pi/4, otherwise:}`. 若干次迭代后, 得到一个近似对角阵 `bm Lambda = bm T^-1 bm A bm T`, `bm Lambda` 的对角元即为 `bm A` 的特征值. 而 `bm T = bm T_1 bm T_2 cdots bm T_n` 的每一列都是 `bm A` 的一个特征向量. 由于 `bm T` 为正交矩阵, 这些特征向量恰好是单位向量. 
# numpy 求特征向量
np.linalg.eig(A)

输出 eigvecs 的每一列表示一个特征向量.

    QR 迭代法 适用于求解一般实方阵 `bm A` 的特征值; 如果 `bm A` 是对称的可逆矩阵, 那么还能求出特征向量. 步骤如下:
  1. 先记 `bm A_1 = bm A`, 将它分解为正交矩阵乘以上三角阵: `bm A_1 = bm Q_1 bm R_1`.
  2. 然后记 `bm A_2 = bm R_1 bm Q_1 = bm Q_1^-1 bm A bm Q_1`, 于是 `bm A_2` 与 `bm A_1` 相似.
  3. 反复迭代使 `bm A_n` 收敛于上三角阵, 此时主对角元为 `bm A` 的特征值.
  4. 特别当 `bm A` 对称时, `bm A_n` 收敛于对角阵, 此时 `bm Q = bm Q_1 cdots bm Q_(n-1)` 的每一列是 `bm A` 的特征向量.
  5. `bm A` 不对称或不可逆时, 上述方法无法求出特征向量. 此时可以将特征值代入到特征方程 `(bm A - lambda bm I) bm x = bb 0` 来求解特征向量.

输出 eigvecs 的每一列表示一个特征向量, 重数大于 1 的特征值可能对应相同的特征向量 (见例3).