Numeric Analysis for Beginners

Chapter 1. Basic Concepts

• 相对误差与绝对误差；

• 求根问题：$\text{For }f:\mathbf{R}\rightarrow\mathbf{R},\text{ find }x^*\text{ such that }f(x^*)=0$；

  假设有估计解 $x_{est}$，但是 $0\lt |f(x_{est})|\ll1$，那么我们也许不知道 $|x_{est}-x_0|$，但是我们一定知道 $|f(x_{est})-f(x_0)|\equiv|f(x_{est})|$；

  • 前向误差：估计解与实际解的差值（就是上面的 $|x_{est}-x_0|$，一般我们不知道）；

  • 后向误差：使得估计值正确所要让 problem statement 改变的 delta（就是上面的 $|f(x_{est})-f(x_0)|\equiv|f(x_{est})|$，一般我们能算出来）；

  • Well-Conditioned（insensitive）：$\text{Small backward error}\Rightarrow\text{Small forward error}$；

  • Poor-Conditioned（sensitive/stiff）：$\text{Small backward error}\nRightarrow\text{Small forward error}$；

  • Condition Number：$CN=\dfrac{\text{Forword Error}}{\text{Backward Error}}$；

    在寻根问题中，很容易得到 $CN=\dfrac{1}{|f^\prime(x^*)|}$；

• Carefully Implementation

  • 防止溢出的方法：Element Scaling（例如求 $||\vec{x}||_2$）；

Chapter 2. Linear System and LU

2.1 Review

• Over Determined / Under Determined / Completely Determined.

  “高” 的矩阵可能是 over determined（更多限制条件）；

  “宽” 的矩阵可能是 under determined（解更可能有无限个）；

• 线性方程组的结论复习：一个线性方程组 $A\vec{x}=\vec{b}$ 有两个不同解 $\vec{x_0},\vec{x_1}$，则它有无穷多解；

  $\vec{x_0},\vec{x_1}$ 都是线性方程组的解，则它们的线性组合都是该线性方程组的解；

• 经验结论：解线性方程组 $A\vec{x}=\vec{b}$ 能不解 $A^{-1}$ 就不解它（计算量和准确性）；

• 初等行变换：左乘初等矩阵；

  初等矩阵是指由单位矩阵经过一次基本行/列变换得到的矩阵。基本行/列变换包括以下三种操作：

  • 交换型初等矩阵：形式为单位矩阵的两行（或列）交换得到的矩阵。例如，在3阶方阵中，交换第一行和第二行的单位矩阵为：

    $$\left[\begin{matrix}0&1&0\\1&0&0\\0&0&1\end{matrix}\right]$$

  • 倍加型初等矩阵：单位矩阵的某一行（或列）乘以非零常数后加到另一行（或列）得到的矩阵。例如，在3阶方阵中，将第一行的两倍加到第三行的单位矩阵为：

    $$\left[\begin{matrix}1&0&0\\0&1&0\\2&0&1\end{matrix}\right]$$

  • 倍乘型初等矩阵：单位矩阵的某一行（或列）乘以非零常数得到的矩阵。例如，在3阶方阵中，将第三行的元素都乘以2的单位矩阵为： $$\left[\begin{matrix}1&0&0\\0&1&0\\0&0&2\end{matrix}\right]$$

  将一个矩阵左乘对应的初等矩阵就是在对行进行对应的变换；

  将一个矩阵右乘对应的初等矩阵就是在对列进行对应的变换；

• 高斯消元：先 forward substitution（将第 i 行的前 i - 1 个元素清零、第 i 个元素置 1），再 back substitution（向之前代入）；

2.2 LU Factorization

LU 分解：当参数矩阵不变，只有 $\vec{b}$ 不同时，我们可以节省重复计算的步骤。这种方法就是 $LU$ 分解。

• 我们先将 $A$ 矩阵使用高斯消元为 $A=LU$（$U$ 是高斯消元的 forward substitution 的结果）；
• 这样可以分开计算 $L\vec{y}=b$ 中的 $\vec{y}$、$U\vec{x}=\vec{y}$ 解得 $\vec{x}$，每一步都是 $O(n^2)$；并且得到上/下三角矩阵的逆要比一般的 $A$ 简单；

2.3 Linear System

• 线性回归预测

• 线性系统拟合曲线：n 次试验构建一个线性方程组；

  

  线性系统拟合非线性曲线：可以让 $f$ 是非线性的，构建一个非线性函数 $f_{ij}$ 的矩阵；

  

  一个特殊的例子是范特蒙德系统 $f=a_0+a_1x+a_2x^2+\cdots$；

  我们还可以借助傅里叶展开 $f=a\cos(x+\phi)$；

• $A\vec{x}=\vec{b}$ 无解情况如何接近？凸优化：$\min\limits_{\vec{x}}||A\vec{x}-\vec{b}||_2^2$；

  

• Tikhonov regularization：这种正则化可以向 under determined 的情况加入限制，有助于防止过拟合、数据抖动、减轻数据噪声影响（对高斯白噪声效果好）等问题：

  $$\min\limits_{\vec{x}}||A\vec{x}-\vec{b}||_2^2+\alpha||\vec{x}||_2^2\quad(0\lt\alpha\le1)\\ \Longrightarrow0=2A^TA\vec{x}-2A^T\vec{b}+2\alpha\vec{x}\\ \Longrightarrow(A^TA+\alpha I_{n\times n})\vec{x}=A^T\vec{b}$$

• 稀疏矩阵存储：变换为低维数据（压缩信息的相关性）；

  • 普通有规律的矩阵，可以通过一些变换转换为稀疏矩阵；

2.5 Cholesky Factorization

注意到重要的矩阵 $A^TA$ 有如下性质：

• 正定对称阵（Positive definite Matrix）其中 $A\ne0$；

而我们又发现正定对称矩阵的 LU 分解非常特殊，$U=L^T$，所以所有正定对称阵可以分解为 $A=LL^T$，这就是 Cholesky 分解；

• 在阶数较小的情况，直接可以使用待定系数法求解 $L$；

• 在阶数较大的情况，使用迭代法求解 $L$，算法如下：

  1. 初始化 $L=0$（全零矩阵，同时是下三角矩阵）；

  2. 对矩阵的每一列 $j$：

     • 先计算该列上的对角元 $L_{jj}=\sqrt{A_{jj}-\sum\limits_{k=1}^{j-1}L_{jk}^2}$（$j\gt1$），其中 $L_{11}=\sqrt{A_{11}}$；

       不难发现，$L$ 每行对角元只与 $A_{jj}$ 和当前行排在 $L_{jj}$ 之前的元素平方和有关；

       可以将 $A$ 分块，即可推出这个结论；

     • 在第 $j$ 列中，继续对第 $j$ 行之后的每一行 $i\gt j$，计算 $L_{ij}=\dfrac{1}{L_{jj}}(A_{ij}-\sum\limits_{k=1}^{j-1}L_{ik}L_{jk})$，

       结论，$L_{ij}$ 与 $A_{ij}$ 和：排在 $L_{jj}$ 前面的元素向量 与 排在 $L_{ij}$ 前面的元素向量的点积有关；

     完成上面两个步骤后，矩阵的第 $j$ 列全部计算完成；

Chapter 3. Norms, Sensitivity & Conditioning in Matrix

3.1 Definitions of Norms in Matrix

引入：在浮点数计算时，如果在处理 $||A\vec{x_0}-\vec{b}||$ 时，它距离 0 有多接近才能相信 $x_0$ 是解？

也就是说，如何衡量 $(A+\delta A)\vec{x}=\vec{b}+\delta\vec{b}$ 求解下 $\vec{x}$ 解的变换幅度？

我们再次引入向量的范数：$||\vec{x}||_p=(\sum\limits_{i=1}^nx_i^p)^{1/p}$；

• 注意到 $p\rightarrow\infty$ 是 $||\vec{x}||=\max\{|x_1|,|x_2|,\ldots,|x_n|\}$；
• $||\vec{x}||=0\quad\text{iff}\quad \vec{x}=0$；
• $||c\vec{x}||=|c|||\vec{x}||,\space c\in\mathbf{R},\vec{x}\in\mathbf{R^n}$；
• $||\vec{x}+\vec{y}||\le||\vec{x}||+||\vec{y}||,\space\forall\vec{x},\vec{y}\in\mathbf{R}$；

我们定义两个范数等价（$||\cdot||_p\equiv||\cdot||_q$） 当且仅当 对于任意 $\vec{x}\in\mathbf{R^n}$，都存在 $c_{low}||\vec{x}||\le||\vec{x}||\le c_{high}||\vec{x}||$（同阶）；

推论：$\mathbf{R^n}$ 上的任意范式等价；

我们再定义矩阵的范数：

• 定义方法 1，“unrolled construction”（元素形式范数，entrywise norm）：将矩阵 $A_{m\times n}$ 按列展开（第 $n+1$ 列排在第 $n$ 列下方），得到一个 $\vec{a}\in\mathbf{R^{mn}}$ 的向量，而向量 $\vec{a}$ 的范数就是 $A$ 的范数；

  $$||A_{m\times n}||_p=(\sum\limits_{i=1}^m\sum\limits_{j=1}^n|a_{ij}|^p)^{1/p}$$
  • 二维元素形式范数又称 “Frobenius Norm”，记作 $||A||_{Fro}$；注意到 $||A||_F=\sqrt{\text{tr}(AA^T)}$；

• 定义方法 2，“induced construction”（诱导范数，又称算子范数，operator norm）：描述了矩阵代表的线性变换 $A$ 对 $\vec{x}$ 作用最长伸展的比例。即：

  $$||A||=\max_{\vec{x}\ne0}\dfrac{||A\vec{x}||}{||\vec{x}||}=\max\limits_{||\vec{x}||=1}||A\vec{x}||$$

  诱导范数的常用结论如下：

  • $||A||_1=\max\limits_{1\le j\le n}\sum\limits_{i=1}^m|a_{ij}|$（1-诱导范数就是一列中的元素模之和，再取最大值）；

  • $||A||_\infty=\max\limits_{1\le i\le m}\sum\limits_{j=1}^n|a_{ij}|$（$\infty$-诱导范数就是一行中的元素模之和，再取最大值）；

  • $||A||_2=\sqrt{\max\limits_{k}\lambda_k}$（2-诱导范数，又称谱范数，是 $A$ 的最大奇异值的开根号。也就是说，当 $A$ 不可逆时就是 $A^TA$ 的最大特征值开根号）；

    可以从图形方法理解：

    

• 定义方法 3，“eigenvalue construction”（schatten 范数，使用矩阵奇异值定义），具体定义较为复杂，不进一步了解；

3.2 Definition of Condition Number in Matrix

我们还要定义一个矩阵的条件数。一般条件数想要看的是一个矩阵构成的线性方程组随误差的变化情况。于是我们需要构建建模一个公式：$(A+\varepsilon\delta A)\vec{x}(\varepsilon)=\vec{b}+\varepsilon\delta\vec{b}$，于是可以得到在误差 $\varepsilon$ 下，

• $\vec{x}(\varepsilon)$ 表示在误差 $\varepsilon$ 下的 $\vec{x}$ 测量值；

• $\dfrac{d\vec{x}}{d\varepsilon}|_{\varepsilon=0}=A^{-1}(\delta\vec{b}-\delta A\vec{x}(0))$ 表示 $\vec{x}$ 随误差的变化率；

• $||\vec{x}(\varepsilon)-\vec{x}(0)||$ 表示前向误差，$\dfrac{||\vec{x}(\varepsilon)-\vec{x}(0)||}{||\vec{x}(0)||}\le|\varepsilon|||A^{-1}||||A||(\dfrac{||\delta\vec{b}||}{||\vec{b}||}+\dfrac{||\delta A||}{||A||})+O(\varepsilon^2)$ 表示归一化的前向误差；

  泰勒展开证明上式；

• 于是我们定义方阵 $A\in\mathbf{R^{n\times n}}$ 的条件数为：$\text{cond}\space A\equiv\kappa\equiv||A||||A^{-1}||$；

  如果 $A$ 不可逆，则条件数为 $\infty$；

  得出结论：$\dfrac{||\vec{x}(\varepsilon)-\vec{x}(0)||}{||\vec{x}(0)||}\le\varepsilon\cdot D\cdot\kappa+O(\varepsilon^2)$；

  可以知道，一个矩阵的条件数描述的性质和行列式不一样，条件数与常系数缩放无关；

• 重要推论 1：$\text{cond}\space A=\dfrac{\max\limits_{\vec{x}\ne0}\frac{||A\vec{x}||}{||\vec{x}||}}{\min\limits_{\vec{y}\ne0}\frac{||A\vec{y}||}{||\vec{y}||}}$（可以使用诱导范数直接推得）；

几何上这么理解：在 $A$ 代表的线性变换下，对任意 $\vec{x}$ 作用伸长最长的比例 与 伸长最短的比例 之比；

由这个几何关系，我们可以推出第二个推论：

• 重要推论 2：$||A^{-1}\vec{x}||\le||A^{-1}||||\vec{x}||$，因为 $||A^{-1}||$ 就是所有向量拉伸最长的比例了；

注：

• 可以知道，一个矩阵的条件数越大，它所构成的线性方程组代表的线性系统对微小扰动越敏感的（解周围小范围变化自变量，总体值变化很大）；而一个矩阵的条件数越接近 1，则这个线性系统对微小扰动越不敏感；
• 由于很难求一个矩阵的逆，因此一般对条件数的讨论是讨论其上下界；

Chapter 4. Column Spaces & QR

考虑特殊矩阵的条件数 $\text{cond}\space A^TA\approx(\text{cond}\space A)^2$（在 $A$ 可逆的情况下）；

$$\begin{aligned} \text{cond}\space A^TA&=||A^TA||\space||(A^TA)^{-1}||\\ &\approx||A^T||\space||A||\space||A^{-1}||\space||(A^T)^{-1}||\\ &=||A||^2||A^{-1}||^2\\ &=(\text{cond}\space A)^2 \end{aligned}$$

所以，我们对于一般的矩阵可以认为 $A^TA$ 越接近单位矩阵 $I_{n\times n}$，$A\vec{x}=\vec{b}$ 更容易解；

此外，$A^TA$ 的计算可以这么理解：

所以等价于我们希望 $A$ 是正交矩阵（各个列向量正交归一），而且恰好正交矩阵代表的变换不改变向量的长度（所以正交矩阵的条件数是 1）；

现在再回来看 $A^TA\vec{x}=A^T\vec{b}$，我们知道这个线性方程组的解就是 $\min\limits_{\vec{x}}||A\vec{x}-\vec{b}||$ 的解，就相当于将 $A$ 拆成列向量 $(\alpha_1,\alpha_2,\ldots,\alpha_n)$，将 $x_1\alpha_1+x_2\alpha_2+\cdots+x_n\alpha_n$ 逼近 $\vec{b}$（将 $\vec{b}$ 投影到 $A$ 列向量做基向量的线性空间上）；

回忆线性代数的性质：

对任何实矩阵 $A\in\mathbf{R^{m\times n}}$ 和可逆方阵 $B\in\mathbf{R^{n\times n}}$，有 $\text{col}\space A=\text{col}\space AB$，即对任意矩阵进行初等行变换不影响矩阵的列空间；

那么有没有办法对 $A$ 一直进行初等行变换，使得 $A$ 变成正交阵？这样 $A$ 的列空间不变（即原问题的解不变），但是 $A$ 成正交阵后非常好求解。

这种方法就是 QR 分解。我们将一般矩阵分解为一个正交阵（$A$ 的一组正交基）和上三角矩阵的乘积（上三角矩阵 $R$ 可以理解为一系列初等行变换）。

显然，对于任意一个 $A$，若 $r(A_{m\times n})=n$（$m\ge n$），则 $A$ 都能进行 QR 分解。

一旦我们将 $A^TA\vec{x}=A^T\vec{b}$ 进行 QR 分解：$A=QR$，那么 $\vec{x}=R^{-1}Q^T\vec{b}$，我们发现 $R$ 上三角矩阵容易求逆，就不需要计算 $A^TA$ 的逆了。

现在，QR 分解有 2 种方法。

• 一种是通过 施密特正交化。这很好理解：

  现在回忆线性代数的 施密特正交化。这就是得到 $A=QR$ 的一种方法：

  1. 先对 $A$ 施密特正交化，再归一化就能得到基向量相同的正交矩阵 $Q$；
  2. 反解出 $R$：因为正交矩阵 $Q^TQ=I$，因此 $R=Q^TA$；

  那么怎么 Gram-Schmidt 正交化？

  对 $A$ 拆成的一组基 $(\alpha_1,\alpha_2,\ldots,\alpha_n)$，可以这么取正交基：

  • $v_1=\alpha_1$，第一个向量随便取；
  • $v_2=\alpha_2-\dfrac{\alpha_2\cdot v_1}{v_1\cdot v_1}v_1$，第二个向量取 $\alpha_2$ 的时候，需要剔除第一个取得的基向量相关方向的分量：$\dfrac{\alpha_2\cdot v_1}{v_1\cdot v_1}$ 就是 $\alpha_2$ 在 $v_1$ 上的投影长度！
  • $v_3=\alpha_3-\dfrac{\alpha_3\cdot v_2}{v_2\cdot v_2}v_2-\dfrac{\alpha_3\cdot v_1}{v_1\cdot v_1}v_1$，第三个向量就减去 $\alpha_3$ 在 $v_1,v_2$ 上的投影分量就行。
  • ……（依此类推）

  最后别忘了归一化。

  不难发现，这种操作实际上在不断地对 $A$ 右乘上三角矩阵 $R_i$（进行初等列变换），使得：$AR_1R_2\cdots R_n=Q$；

  相对地，下面的 Householder 变换方法，就是不断地进行正交变换（特别地，镜像变换），调整某一列的其他元素为 0，即不断地对 $A$ 左乘 Householder 矩阵（一种正交阵）$H_i$，使得：$H_nH_{n-1}\cdots H_1=R$；

• 另一种是通过 Householder 变换（这是一个著名的变换，它代表了镜像变换，显然是一种正交变换），它所对应的矩阵就是 Householder 矩阵。

  如图：

  

  假设已知一个向量 $\eta$，想要关于某个法线方向 $l$ 对称。为了方便，我们记与 $l$ 正交的一个从 $\eta$ 一边指向另一边的单位向量为 $\omega$，则由几何关系可知对称后的向量 $\xi$ 满足：$\xi-\eta=2\omega(\omega^T\xi)$，其中 $\omega^T\xi$ 就是 $\omega$ 和 $\xi$ 的点积（投影长度）；

  因此 $\eta=(I-2\omega\omega^T)\xi$，我们发现，这镜面变换的矩阵就是 $H=I-2\omega\omega^T$，左乘它会将列向量变换到 $\omega$ 对应法线的 $\omega$ 指向的另一侧。

  紧接着，我们发现这个矩阵 $H$ 有这个性质：

  若 $H$ 为 Householder 矩阵，则 $\left[\begin{matrix}I_r&0\\0&H\end{matrix}\right]$ 也是 Householder 矩阵；

  于是我们可以这么进行迭代：

  • 写 $A$ 的基向量 $(\alpha_1,\alpha_2,\ldots,\alpha_n)$，做第一次 householder 变换，使得 $\omega_1=\dfrac{\alpha_1-||\alpha_1||_2\cdot e_1}{||\alpha_1-||\alpha_1||_2\cdot e_1||_2}$，得到第一个 householder 矩阵 $H_1=I-2\omega_1\omega_1^T$，这样 $H_1A$ 的第一列除了第一个元素全部归 0：

    $\omega_1$ 可以理解为 $\alpha_1$ 减去在 $e_1$ 的分量（要镜像的方向已经得到了），再归一化，得到单位的镜像向量；

    $H_1\alpha_1$ 的变换就将 $\alpha_1$ 变换到与 $e_1$ 同一个方向上了。

    

  • 记 $H_1A$ （注意是变换后的矩阵）关于 $a_{11}$ 的余子式为 $B_1$，则写 $B$ 的基向量 $(\beta_1,\beta_2,\ldots,\beta_{n-1})$，那么同理 $\omega_2=\dfrac{\beta_2-||\beta_2||_2\cdot e_1}{||\beta_2-||\beta_2||_2\cdot e_1||_2}$，得到第二个 householder 矩阵 $H_2=I-2\omega_2\omega_2^T$，这样 $H_2H_1A$ 的第一列和 $H_1A$ 一样、第二列除了前两个元素，后面的元素全部归 0:

    

    注意，$e_1$ 总是第一个元素为 1、其他元素为 0 的单位向量；

  • 重复上面的操作，最后 $H_{n-1}H_{n-2}\cdots H_1A=R$，$R$ 是一个上三角矩阵：

  • 这样，我们得到了 $A$ 分解出的 $R$，最后反代 $Q=H_1H_2\cdots H_{n-1}$（注意正交阵的性质）；

Chapter 5. Eigenvalue & Eigenvector

5.1 Overview

在采集多维数据时，需要考虑各个维度间的相关性，以降低数据的维度。

举个例子，假设有组（该组有 $m$ 个数据） $n$ 维数据 $(v_1,v_2,\ldots,v_m),\space v_i\in\mathbf{R^n}$。

如果我们只知道某些维度上的确切数据，于是我们就像想将任意一个 $n$ 维数据用某几个维度去拟合整体数据。这样可以非常方便地讨论数据的整体特性。

这样，数据矩阵的特征值就能派上用场了。

除了这个问题，还有其他一些问题可以借助特征值进行解决，例如：

• Optimize $||A\vec{x}||_2$，固定 $||\vec{x}||_2=1$；
• ODE/PDE（常微分方程、偏微分方程）的近似解：$\vec{y}^\prime=B\vec{y}$；
• Rayleigh quotient（瑞利商）：$\dfrac{\vec{x}^TA\vec{x}}{||\vec{x}||_2^2}$；

回忆下线性代数中对于特征向量/特征值的重要结论：

• 每个 $n$ 阶方阵至少有一个特征向量（复向量），最多有 $n$ 个不同的特征值；
• 对应不同特征值的特征向量是线性无关的；

5.2 Review: Diagonalizable Matrix

这里再复习一下矩阵对角化的知识：

矩阵对角化的意义？

• 可快速计算 $A^k$；

• 可计算 Markov 过程中的平稳分布 $\pi$；

• 可计算差分方程 $u_{k+1}=Au_k$ 描述的离散动力系统的长期行为；
• ……

矩阵对角化的方法？

1. 求出矩阵 $A$ 的所有特征值 $\lambda_i$；

2. 通过 $A$ 的每个特征值，以及特征值的代数重数，来判断 $A$ 是否可对角化。具体来说：

   代数重数就是在判断特征值重复的次数、几何重数就是在描述特征向量重复的维数（就是零空间的维数）。注意每个特征值的几何重数一定小于等于代数重数（因为对应不同特征值的特征向量是线性无关的，而特征值可以重复）。

   这里 $A$ 要可对角化，就必须满足下面两种情况之一：

   • $A$ 的所有 $n$ 个特征值互不相等（代数重数 $n$）。而由于对于不同特征值的特征向量必然线性无关，所以几何重数一定也为 $n$；
   • $A$ 所有重根下，$k$ 重特征值是否有 $k$ 个线性无关的特征向量。这里就是在要求这个代数重数为 $k$ 的特征值的几何重数是不是也是 $k$；

   所以上面的两个要求总体在说：$A$ 的几何重数和代数重数是否相等？

   如果相等，表示 $A$ 所代表的线性变换没有改变被变换方的维度，因此 $A$ 可以分解为 $n$ 个线性无关的正交基向量。也就是可以相似对角化。如果不满足，则无法相似对角化。

3. 最后将特征向量与特征值对应起来：$\Lambda=\text{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n)$，$P=(\alpha_1,\alpha_2,\ldots,\alpha_n)$，则 $P^{-1}AP=\Lambda$；

5.3 Definition of Spectrum Radius in Matrix

引入一个新的定义：矩阵的谱半径（或称 “矩阵的谱”）。

谱半径的意义？

• 估计一个矩阵的特征值；
• 计算一个不可逆矩阵的广义逆矩阵；
• ……

谱半径的计算？

$\rho(A)=\max\limits_i|\lambda_i|$（矩阵 $A$ 的谱半径等于其所有特征值的模的最大值。注意特征值包含复数！）；

谱半径和范数的关系？

谱半径和矩阵范数一样，都是矩阵的函数：$f:R^{m\times n}\rightarrow R$；

但是它们二者本质上真的不一样，一定要和 2-诱导范数（也就是谱范数）区分开。

二者间有一些重要结论：

• 任意复数域上的矩阵 $A$，其谱半径 $\rho(A)$ 不大于 $A$ 的任何一种诱导范数，即：$\rho(A)\le||A||$；

  含义：矩阵的谱半径是其任意一种范数的下界；

  作用：使用方便求解的范数对谱半径进行估算；

• Gelfand 定理：$\rho(A)=\lim\limits_{k\rightarrow\infty}||A^k||^{1/k}$；

  • 矩阵序列 $I,A,A^2,\ldots,A^k,\ldots$ 收敛于 0 的充要条件：$\rho(A)\lt1$；
  • 级数 $I+A+A^2+\cdots$ 收敛于 $(I-A)^{-1}$ 的充要条件：$\rho(A)\lt1$；

5.4 Extend to $\mathbf{C^{m\times n}}$

现在将线性空间扩展到复数域，我们多出如下定义：

• 共轭转置 $A^H$；

• 厄密矩阵（Hermitian Matrix）：$A^H=A$；

  注意和实对称矩阵不一样。

  由量子力学的厄密算符可以得到如下所有结论（量子力学考题）：

  • 厄密矩阵所有特征值为实数；
  • 厄密矩阵属于不同特征值的特征向量相互正交；
  • ……

  厄密矩阵因为是复数域上，和实数域上的实对称矩阵很相似（不如说实对称矩阵是厄密矩阵的特殊情况），所以厄密矩阵和实对称矩阵一样，二者一定可以相似对角化（几何重数一定等于代数重数，或者说一定有 $n$ 个线性无关的基向量）；

  它们相似对角化很简单：$A=X^{-1}\Lambda X=X^T\Lambda X$（显然 $X$ 是正交矩阵）；

再引入一些 “奇怪” 的运算：

对于一个 半正定/正定 的对称矩阵 $A\in S_{+}$，定义其平方根：$A^{1/2}$，因为一定能找到 $P$ 使得 $P^2=A$；

5.5 Application: Use Matrices to Solve Problems

• 例如 $\vec{y}^\prime=\lambda\vec{y}$，可以看成一个求导变换 $D\vec{y}=\lambda\vec{y}$，求 $D$ 的特征向量就是 $\vec{y}$ 的解；

• 照片曝光的例子：数据集中有 $n$ 个数据，$\omega_{ij}\ge0$ 表示第 $i$ 和 $j$ 数据集之间的某个指标的相似性，$\omega_{ij}=\omega_{ji}$；我们想将这些数据集以这个指标 $x_i$ 衡量起来，要求相似性越高的数据，$x_i$ 的值也应该相近；

  为了完成这个任务，可以定义一个目标函数 $\sum\limits_{ij}\omega_{ij}(x_i-x_j)^2$，对它最小化优化就行。

  但是需要一些限制条件，防止 $x_i\equiv const$ 无意义的情况，例如 2-范数为 1 $||\vec{x}||_2^2=1$（归一化）、$\vec{1}\cdot\vec{x}=0$（指标 $x_i$ 均值为 0，方便统计）；

  所以目标函数可以简化为：$2x^T(A-W)x$，其中 $W=(\omega_{ij})_{n\times n}$；找到 $A-W$ 的第二小特征值（最小特征值是 0，已经被限制条件排除了）对应的特征向量就是解。

• 计算 $A^k$：

对一个实对称阵 $A$，假设其特征值 $\lambda_1,\ldots,\lambda_n$ 从大到小排列（$\lambda_{i+1}\ge\lambda_{i}$），那么它由实对称阵特征向量的完备性，我们可以用 $A$ 特征向量 $(x_1,\ldots,x_n)$ 来表示任意 $n$ 维向量：$A\vec{v}=c_1A\vec{x_1}+\cdots+c_nA\vec{x_n}=c_1\lambda_1\vec{x_1}+\cdots+c_n\lambda_n\vec{x_n}$（和量子力学将力学量使用它对应的算符的本征函数展开是一样的）；

因此我们发现，$A^k$ 对 $\vec{v}$ 作用的效果就取决于最大的特征值及其特征向量了：

$A^k\vec{v}\approx c_1\lambda_1^k\vec{x_1},\space\text{assume that}\space c_1\ne0,\lambda_1\gt\lambda_2$（$k$ 要求足够大）；

问题是，如果 $|\lambda_1|\gt1$ 时，$A^k\vec{v}\rightarrow\infty$，所以每次迭代都做一次归一化即可：

$\vec{\omega_k}=A\vec{v_{k-1}},\space\text{ where }\vec{\omega_{k}}=\dfrac{\vec{v_{k-1}}}{||\vec{v_{k-1}}||}$；

所以通过这个 power iteration 方法我们就能估计出 $A$ 的最大特征值；

又注意到 $A$ 特征值的倒数，正好是 $A^{-1}$ 对应的特征值：$A\vec{v}=\lambda\vec{v}\Rightarrow A^{-1}\vec{v}=\dfrac{1}{\lambda}\vec{v}$；

那么对 $A^{-1}$ 进行 power iteration，就能得到 $A$ 的最小特征值。

$A^{-1}$ 的 power iteration 可以借助 $LU$ 分解加速。

另外，由于正确结果收敛较慢，因此我们可以使用 “shifted inverse iteration”：

$A\vec{v}=\lambda\vec{v}\Rightarrow(A-\sigma I)\vec{v}=(\lambda-\sigma)\vec{v}$，可以得到如下的迭代过程（猜测 $\sigma_k\approx\lambda_k$）：

$$\vec{\omega_k}=(A-\sigma_k I)^{-1}\vec{v_{k-1}},\space\vec{v_k}=\dfrac{\vec{\omega_k}}{||\vec{\omega_k}||},\space \sigma_{k+1}=\dfrac{\vec{v_{k}}^TA\vec{v_k}}{||\vec{v_k}||_2^2}$$

5.6 Similarity Transformations

• 借助 QR 分解进行相似变换：$A=QR$，则 $Q^{-1}AQ=RQ$，所以可以迭代：$A_{k+1}=R_kQ_k,\space A_k=Q_kR_k$；

  好的结论：$A_\infty=Q_\infty R_\infty=R_\infty Q_\infty$；

5.7 SVD (Singular Value Decomposition)

回忆一个矩阵的诱导范数：

$$||A||=\max_{\vec{x}\ne0}\dfrac{||A\vec{x}||}{||\vec{x}||}=\max\limits_{||\vec{x}||=1}||A\vec{x}||$$

注意到 $||\vec{v}||=1$ 时，$||A||^2=||A\vec{v}||^2=\vec{v}A^TA\vec{v}$，因此只需要

特征值只有方阵才能讨论。有没有一种研究矩阵更普遍特征性质的分解呢？它就是奇异值分解。

可以证明，任何矩阵可以分解为：$A=U\Sigma V^{-1}$，其中 $U,V$ 为正交矩阵，$\Sigma$ 为对角矩阵（可以不是方阵）；

在物理上 $U,V$ 表示旋转变换（rotation），$\Sigma$ 表示伸缩变换（scale）；

那么 $U$、$V$ 代表什么？

我们发现：$A^TA=V(\Sigma^T\Sigma)V^T$，显然 $\Sigma^T\Sigma$ 是个对角方阵。而 $A^TA$ 有很好的性质：它是半正定、对称矩阵。

因此我们惊喜地发现，这就相当于对 $A^TA$ 完成了相似对角化。$\Sigma^T\Sigma$ 对角元存放的是 $A^TA$ 的特征值，$V$ 的列存放的是 $A^TA$ 对应的特征向量；

同理，$AA^T=U(\Sigma^T\Sigma)U^T$，所以，我们得出以下结论：

• $U$ 和 $V$ 分别是 $AA^T$ 和 $A^TA$ 的归一化特征向量组成的正定矩阵；
• $U$ 和 $V$ 特征值按序相同（从大到小排列），对于实矩阵而言都大于等于 0；

SVD 分解还可以写成一系列 $U,V$ 向量的外积线性组合：$A=\sum\sigma_i\vec{u_i}\vec{v_i}^T$；

SVD 可以有哪些用处？比如定义一个一般矩阵的 “伪逆”（pseudo inversion）：

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

注意，对不一定为方阵的对角阵的伪逆：$\Sigma^{+}$ 就是将 $\Sigma$ 对角元的元素求倒数，放在对角位置，并且转置矩阵的长宽。

伪逆有一些很好的性质：

• 当 $A$ 为可逆方阵时，$A^+=A^{-1}$；
• 当 $A$ 的行秩大于列秩时（overdetermined），$A^+\vec{b}$ 给出了最小二乘结果；
• 当 $A$ 的列秩大于行秩时（underdetermined），$A^+\vec{b}$ 给出了 $A\vec{x}\approx\vec{b}$ 的最小二乘结果；

5.7.1 Application: Orthogonal Procrustes Theorem

考虑一个问题：将一组向量 $A$ 通过正交变换的方式映射到新的一组向量 $QA$，让这组新的向量与给定的一组向量 $B$ 的差异尽可能的小（使用 Frobenius 范数衡量）。这在图像处理领域比较常用。

用数学方法表达就是，求正交阵：$\hat{Q}=\arg\min\{||QA-B||_{F}\}$（易知，$A,B$ 规模相同）；

由于 Frobenius 范数可以表达为 $||C||_F=\sqrt{\mathrm{tr}(C^TC)}$，故上式可计算：

$$\begin{aligned} \hat{Q}&=\arg\min\{\mathrm{tr}((QA-B)^T(QA-B))\}\\ &=\arg\min\{\mathrm{tr}(((QA)^T-B^T)(QA-B))\}\\ &=\arg\min\{\mathrm{tr}((A^TQ^T-B^T)(QA-B))\}\\ &=\arg\min\{\mathrm{tr}(A^TQ^TQA-BA^TQ^T-B^TQA+B^TB)\}\\ &=\arg\min\{\mathrm{tr}(A^TA+B^TB-A^TQ^TB-(A^TQ^TB)^T)\}\\ &=\arg\min\{\mathrm{tr}(A^TA)+\mathrm{tr}(B^TB)-2\mathrm{tr}(A^TQ^TB)\}\\ \end{aligned}$$

注意到前两项在确定问题时就已知，所以最小化问题直接转换为最大化问题：

$$\begin{aligned} \hat{Q}&=\arg\max\{\mathrm{tr}(A^TQ^TB)\}\\ &=\arg\max\{\mathrm{tr}(Q^TBA^T)\}\\ \end{aligned}$$

（注意 $\mathrm{tr}(PQ)=\mathrm{tr}(QP)$）

下面利用 SVD 分解做一个巧妙变换（假设 $BA^T=U\Sigma V^T$）：

$$\begin{aligned} \hat{Q}&=\arg\max\{\mathrm{tr}(Q^TBA^T)\}\\ &=\arg\max\{\mathrm{tr}(Q^TU\Sigma V^T)\}\\ &=\arg\max\{\mathrm{tr}(V^TQ^TU\Sigma)\}\\ \end{aligned}$$

这样 $V,Q,U$ 都是正交阵，所以 $V^TQ^TU$ 也是正交阵。

可以简单地证明，$\mathrm{tr}(QA)$ 最大（$Q$ 正交阵）时，$Q=I$；

此时 $\hat{Q}^T=VU^T$，$\hat{Q}=UV^T$；

所以这里我们利用 SVD 解出了正交普鲁克定理：

当 $Q=UV^T$（其中 $BA^T=U\Sigma V^T$）时，$||QA-B||_F$ 有最小值。

5.7.2 Application: Principal Component Analysis (PCA)

再考虑一个问题，对于相当大的一个数据集，它包含很多维度，我们基于以下目的需要降维：

• 使得数据集更易使用；
• 降低算法的计算开销；
• 去除噪声；
• 使得结果容易理解；

有一种降维方法就是主成分分析方法（PCA），其主要思想是将 $n$ 维特征映射到 $k$ 维上，这k维是全新的正交特征也被称为主成分，是在原有 $n$ 维特征的基础上重新构造出来的 $k$ 维特征；

可以证明：

最小化 $||X-CC^TX||_F$（其中 $C\in\mathbf{R^{n\times k}}$，$C^TC=I_{k\times k}$，$C$ 是 $U$ 的前 $k$ 列向量，$X=U\Sigma V^T$）取得的 $C$ 就是 $X$ 的主成分。

Chapter 6. Non-linear System

6.1 Root Finding

作出一些假设：

• 连续性；

  连续函数满足的定理：中值定理；

• Lipschitz 特性：绝对值增长速率不快于一阶线性函数；

• $k$ 阶导存在且连续；

方法一：二分法（bisection），利用中值定理锁定根的区间（高中内容），直到根的区间小于一定范围就停止迭代。

• 优点：无条件收敛（unconditionally converge）；

  收敛速度？指数速度减小。

  $|x-x^*|\lt E_k$，其中 $E_k$ 为第 $k$ 轮迭代时的区间宽度（$E_k\le\dfrac{1}{2}E_{k-1},\space E_k=|r_k-l_k|$）；

• 缺点：对函数性质要求严格。

方法二：不动点法（fixed point），通过迭代求解 $g(x^)=x^$ 来得到 $f(x)=g(x)-x$ 的零点。

怎么迭代？方法比较多，但是常用的是最简单的策略：

simple strategy：将 $g(x_{k-1})$ 作为下一轮迭代的 $x_k$ 的值，直至 $|g(x_k)-x_k|$ 小于一定范围；

• 优点：计算简单；

• 缺点：$g$ 必须满足 Lipschitz 特性，或者在根 $x^*$ 及迭代范围附近满足 Lipschitz 局部特性，否则迭代发散。

  对于满足 Lipschitz 特性的情况，$E_k\le cE_{k-1}$（linear）；

  对于其他一般情况：$E_k=|x_k-x^|=|g(x_{k-1})-g(x^)|\le\dfrac{1}{2}(|g^{\prime\prime}(x^*)|+\varepsilon)E_{k-1}^2$（quadratic）；

方法三：牛顿法（Newton’s method），这个方法作出了一个假设，认为函数在零点附近近似线性，可以给出一个猜测值 $x_0$，求该点处切线 $l_0$，取得 $l_0$ 与 x 轴交点为 “更接近零点的点” $x_1$，重复迭代直至稳定。

$$x_{k+1}=x_k-\dfrac{f(x_k)}{f^\prime(x_k)}$$

可以看作求 $g(x)=x-\dfrac{f(x_k)}{f^\prime(x_k)}$ 的不动点。

• 优点：在简单情况下收敛很快；
• 缺点：
  1. $g$ 需要满足局部 Lipschitz 特性，否则不收敛（收敛速度和不动点法同理）；
  2. $f^\prime(x^*)\ne0$，否则永远无法得到正确解；
  3. 某些函数难以求导数。

方法四：割线法（secant method），利用两相近点间的割线近似为切线的思想（$f^\prime(x)\approx\dfrac{f(x_1)-f(x_2)}{x_1-x_2},\space |x_1-x_2|\rightarrow0$），可以借助两次猜测的点的连线（割线）视作切线：$f^\prime(x_k)\approx\dfrac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}$，结合牛顿法解决：

$$x_{k+1}=x_k-\dfrac{f(x_k)}{f^\prime(x_k)}\approx x_k-\dfrac{f(x_k)(x_k-x_{k-1})}{f(x_k)-f(x_{k-1})}$$

• 优点：计算稍微简单一点；
• 缺点：和牛顿法一样存在收敛问题。