数值分析

Introduction

使用计算机计算和现实计算有一个显著差别：计算机计算的精度是有限制的。我们有不同的数据结构，带来不同的计算精度。数值分析这门课要求我们求出精度足够好的结果。它会告诉我们一些近似算法，同时也告诉我们，什么时候它们能用，什么时候又不能用。

我们轻而易举地就能想清楚加减乘除的原理。然而，我们可能没有注意， $\sin, \cos, \tan, \ln$ 是怎么实现的？数值分析这门课可以告诉我们背后的原理。

一个自然的想法：泰勒展开

数学上的准备

怎么求一个积分？简单的想法是把被积函数泰勒展开转化成多项式，然后利用其便于求积分的特性求解。然而，泰勒展开往往是无限的，这就需要我们规定一个提前结束的时机，这也引入了一个相当于泰勒余项(Remainder)的误差。这个误差就是我们需要注意的。

而且，一个棘手的点是：我们没有办法直接和真值进行比较！

另一个棘手的点是：每一个数值本身都与它自己的真值有误差，这也引入了新的误差！我们需要为每一个数据保留多少位，这也相应地成为了一个问题。

误差 Errors

• Truncation Error：与时间有关的一个误差。它代表近似数学引入的误差。

显式的操作，经典的就是一个 for 循环

• Roundoff Error：与空间有关的一个误差。它代表数字在计算机中的表达和数字本身的误差。

数据背后的误差

为了讨论方便，我们用十进制来表示数字。对于一个不能直接以有限数位表示的数字，有两种方式处理精度：

• 四舍五入(Rounding)
• 直接砍掉精度数位后面的位数(Chopping)
• 往上取

大学老师给分(不是)

误差也分为绝对误差(absolute error)和相对误差(relative error)

而有效数字(siginificant digits)应该是一个相对误差概念。可以注意到，对应数值的部分被移到了科学计数法中的指数部分。

有了有效数字的概念，我们应该把数字 0.123 看成 $0.123\pm\epsilon_1$ 使用四舍五 入的方法时，有效数字 0.1 的相对误差是 50%

那么，有一个现象，就是两个相近的数字相减之后，有效位数减小时，相对误差会显著增加。

将一个数除以一个很小的数，绝对误差会放大。这一点是符合直觉的。一个想法是把绝对误差看成两个数字的函数，然后取关于分母的导数，如果导数很大，那么就说明分母的轻微变化会带来较大的变化，也就是较大的误差。对于相对误差，也可以尝试以这样的方法进行分析。

计算的约化也要注意，计算机是每一个单元运算都会进行约化的，不能直接对最终结果进行约化。

也正是因此，虽然对单个数字而言，Rounding 会更精确，但是对于一个一连串的算式而言，未必。

我们此时也自然地想到，解不是确定值，也可能不是个区间，而应该是一个概率密度分布中的某个可能值。

也还是因此，对于同一个算式的不同表达，比如将一个算式通过一些方式结合，分配，也可以导致不同的误差。减少乘除的次数，有可能减少误差。减少单元计算的数量，也有可能减小误差。

数学上等价，不等于数值分析方法上等价！

我们可以手动求导来分析误差，也可以用一些自动求导的工具来分析计算式的误差。

算法和收敛 Algorithms and Convergence

当一个算法中，原始数据的较小变化只引起较小的终解的变化时，它是 stable 的；否则，它是 unstable 的。当它对于某些原始数据 stable 时，它是 conditionally stable 的。

设初始误差为 $E_1$，当做连续的 n 次操作时，如果误差 $E_n$ 约为 $E_1$ 的常数倍时，说误差的增长是 linear 的。如果 $E_n$ 约为 $E_1$ 的以常数为底的指数倍时，则说 误差的增长是 exponent 的。

一元方程的求解

二分法

能使用二分法的前提是可排序。比如说，复数就很难使用二分，因为它没有既定的序。

二分法一定要除以二吗？不一定。只要区间收敛就好。

• 取中间位置时应该用 $p=a+\frac{b-a}{2}$。
• 判断取左边界还是右边界时，需要用 sign()，而非直接相乘看是否小于 0。
• 函数值需要考虑溢出问题。比如说，一个指数函数就容易出现溢出，即使它的横坐标是没溢出的。

如果二分法的区间取的过大，可能会忽略函数的根，如一个先增后减的区间，两侧都小于零，就会误认为中间没解了。

不动点迭代

把方程的根转化成一个等效的等式：$f(x)=0 \Leftrightarrow x=g(x)$。$g$ 的不动点是 $f$ 的根。

不动点定理：Self

令 $g$ 是一个在 $[a, b]$ 连续的函数，且在其中 $g(x)\in[a, b]$。若对于它的导函数 $g'$，存在一个常数 $k\in(0, 1)$ 使得开区间内任意的 $x$，有 $|g'(x)|\le k$，则对任意的 $p_0\in[a, b]$，序列 $p_n=g(p_{n-1})$ 收敛到唯一的不动点 $p\in[a, b]$。

• 利用中间值定理，证明存在不动点。
• 利用中值定理，证明不动点唯一。
• 利用值域被定义域包含的条件，保证迭代过程中，$g(x)$ 的结果始终在定义域中。
• 利用中值定理，证明迭代确实是收敛于那个存在且唯一的不动点。

这里的存在一个 $k\in(0, 1)$ 很重要，这把 $g'(x)$ 和 1 隔开了，避免了极限为 1 的情况。

看证明是有利于记忆定理的。在看证明的过程中，可以理解每一个条件为什么被需要。

数值分析中有一个特点，给出的公式往往是充分的。顶层应用时，即使没有满足定理的条件，但我们还是有可能选择相信它。而在底层处理时，我们才会格外注意它的必要性。

不动点定理：Corollary

如果 $g$ 满足了不动点定理，那么迭代误差的范围为：

$$|p_n-p|\le\frac{1}{1-k}|p_{n+1}-p_n|\quad, \quad |p_n-p|\le\frac{k^n}{1-k}|p_1-p_0|$$

牛顿法

牛顿法也属于一种不动点方法。它的思路是把一个非线性函数线性化。

定理

若 $f$ 在 $[a, b]$ 二阶连续，且存在 $p\in[a, b]$ 使得 $f(p)=0, f'(p)\ne 0$，那么 就存在一个 $\delta > 0$，使得任意初值 $p_0\in[p-\delta, p+\delta]$ 都可以满足 $p_{n+1}=p_n-\frac{f(p_n)}{f'(p_n)}$ 收敛到 $p$。

• 由 $f'(p)\ne0$ 知 $g(x)=x-\frac{f(x)}{f'(x)}$ 在 $p$ 的邻域中连续。

• $g'(x)=\frac{f(x)f''(x)}{f'(x)^2}$ 在 $p$ 的邻域中趋于零，只要 $f''(x)$ 连续 且有限。

• 利用不动点定理知存在一个邻域使得收敛。

作业-2

P54-T13

找到一个迭代次数的范围，使得使用二分法解 $[1, 2]$ 上的方程 $x^3 -x-1=0$ 时，解 的精确度有 $10^{-4}$。同时，给出这个解。

此题主要考虑的是二分法的误差范围。当对 $[a, b]$ 进行二分时，第一次二分迭代的结果 为 $a+\frac{b-a}{2}$，$\epsilon\le\frac{b-a}{2}$。易知，第 $n$ 次二分迭代的误差 范围为 $\epsilon\le\frac{b-a}{2^n}$，此处只需取 $n=14$。具体求解过程略。

P54-T15

令 $p_n=\sum_{k=1}^{n}\frac{1}{k}$，证明 $p_n$ 是发散 diverge 的，即使 $\lim_{n\rightarrow\infty}(p_n-p_{n-1})=0$。

这算是很经典的题目了，典型的方法有积分放缩。这里则采取证明连续多项的和拥有固定下界的方法：

因此总能找到 $p_{2k-1}>p_{k}+\frac{1}{2}$，故 $p_n$ 是发散的。

P64-T3

它们被用于计算 $21^{\frac{1}{3}}$。根据收敛速度，将它们排序。假设 $p_0=1$。

a. 这是一个比较普通的不动点方法，直接对

求导，得

这看起来有些危险，我们不妨先迭代一次，得到 $p_1=\frac{41}{21}$，如果把这看成新的迭代起点，那么可以发现 $f'(x)$ 的范围被控制住了，且拥有明确的上下界。那么，根据不动点定理的推论，可以得出误差的范围为

b. 通过观察可以发现，这是牛顿法，原函数即为

关于牛顿法，有这样一件事情：

牛顿法的迭代公式为

$$p_{n+1}=p_n-\frac{f(p_n)}{f'(p_n)}$$

设迭代终点为 $r$，记 $\epsilon_n=p_n-r$，则根据泰勒展开式，我们有：

$$\epsilon_{n+1}=\epsilon_n-\frac{f'(r)\epsilon_n+\frac{1}{2}f''(r)\epsilon_n^2+O(\epsilon_n^3)}{f'(r)+f''(r)\epsilon_n+O(\epsilon_n^2)}=\frac{f''(r)}{2f'(r)}\epsilon_n^2$$

这也就是所谓的牛顿法具有二次收敛的性质，这个收敛速度已经很快了。

本选项就是一个二次收敛的样例。

需要注意的是，当出现重根， $f'(x)=0$，此时二次收敛不再成立。

c. 看起来很复杂，我们直接尝试代入，然后发现 $p_1=0, p_{k>1}=0$，可知压根收敛不到 $21^{\frac{1}{3}}$。

d. 可以设 $p_n=21^{k_n}$，则有

可知 $k_n$ 的误差为 $O(\frac{1}{2}^n)$，那么

故误差为 $O(\frac{1}{2^n})$。

综上，排序为 $b>d>a>c$。

P65-T19

---

1. 利用不动点定理，证明当 $x_0\gt\sqrt{2}$ 时，

$$x_n=\frac{1}{2}x_{n-1}+\frac{1}{x_{n-1}}, \qquad \text{for } n\ge 1$$

收敛到 $\sqrt{2}$，

2. 利用“当 $x_0\ne\sqrt{2}$，$0< (x_0-\sqrt{2})^2$”，证明：当 $0\le x_0\le\sqrt{2}$，有 $x_1\gt\sqrt{2}$。
3. 利用 (1), (2) 的结果，证明 (a) 中的式子收敛到 $\sqrt{2}$，只要 $x_0\gt 0$。

1. 记 $f(x)=\frac{1}{2}x+\frac{1}{x}$，当 $\sqrt{2}<x\le x_0$ 时，有 $\sqrt{2}<\frac{1}{2}x+\frac{1}{x}<x_0$，即 $f(x)\in(\sqrt{2}, x_0)$。函数连续性显然；再由 $f'(x)=\frac{1}{2}-\frac{1}{x^2}\in(0, \frac{1}{2})$ 知存在 $k\in(0, 1)$。综上，已满足不动点定理。
2. 题意估计是要我们配方，不过基本不等式易得了。
3. 首先，利用 (2) 中结果，当 $x_0\in(0,\sqrt{2})$ 时，不妨先迭代一次，把迭代结果看成新的 $x_0>\sqrt{2}$。然后，由 (1) 易得收敛至 $\sqrt{2}$。 $x_0\in(\sqrt{2}, \infty)$ 时亦然。

对迭代方法的误差分析

假设 $\{p_n\}$ 是一个收敛到 $p$ 的序列，如果存在正常数 $\alpha$ 和 $\lambda$ 使得

$$\lim_{n\rightarrow \infty}\frac{|p_{n+1}-p|}{|p_n-p|^\alpha}=\lambda$$

则说 $\{p_n\}$ 对 $p$ 的收敛阶至少为 $\alpha$

比如说，牛顿法就是一个至少二阶收敛的方法(当不存在重根)。

下面是一个可以确定收敛阶阶数的定理：

设 $p$ 是 $g(x)$ 的一个不动点，若在 $p$ 的邻域内 $g(x)$ $\alpha \ge 2$ 阶连续，且满足

$$g'(p)=g^{(2)}(p)=\ldots g^{(\alpha-1)}(p)=0, g^{(\alpha)}\ne 0$$

则 $g(x)$ 的不动点迭代的收敛阶为 $\alpha$。

重根

对计算机来说，解 $x^2=0$ 远比 $x^2-4x+3=0$ 困难。这是因为，前者要考虑重根，而且对于计算机来说，这实际上是在求 $x^2=\pm \epsilon$，当 $\epsilon<0$，会出现未定义的情况。所以说，解方程有可能是不稳定的。这里的一种处理方法是，改而求解 $(x^2-0)^2=0$。

前面我们提到，对于牛顿法而言，重根的存在会导致它的收敛阶低于 $2$。不过，我们可以做一个修正，也就是构造另一个函数，使得它的根和原函数的根相同，但是在原重根的邻域内不存在重根。

关于邻域内不存在重根这一点，可以通过把函数 $f(x)$ 写成

$$f(x)=(x-p)^k\phi(x)\quad, \quad\phi(p)\ne 0$$

来得出。

Aitken’s Δ² Method

对于一个一阶收敛的迭代序列，有

$$\frac{p_{n+1}-p}{p_{n}-p}\approx \frac{p_{n+2}-p}{p_{n+1}-p}$$

据此，我们可以解出 $p$ 的一个近似解

$$p\approx\frac{p_{n+2}p_n - p_{n+1}^2}{p_{n+2} - 2p_{n+1} + p_n}=p_n-\frac{(p_{n+1}-p_n)^2}{p_{n+2}-2p_{n+1}+p_n}$$

如果引入差分的话，就可以写成

$$p\approx p_n-\frac{(\Delta p_n)^2}{\Delta^2 p_n}$$

于是，我们在使用一阶收敛的方法时，就可以使用 Aitken’s Δ² Method 来加速收敛。

有人可能会考虑把 Aitken’s Δ² Method 加速完的结果再作为新的迭代起点，然而这样破坏了 linear convergence 的条件，是不可取的。

作业-3-1

P86-T11

不动点迭代方法

$$p_{n+1}=g(p_{n})=p_n-\frac{f(p_n)}{f'(p_n)}-\frac{f''(p_n)}{2f'(p_n)}\left[\frac{f(p_n)}{f'(p_n)}\right]^2$$

拥有性质

$$g'(p)=g''(p)=0$$

这会得到立方收敛的结果。比较一下平方收敛和立方收敛。

设另一个平方收敛的序列为 $\{q_n\}$，为了便于分析，我们设它们都收敛到 $0$(其实就是移动坐标轴)，且

可以得出

求解线性方程组 Linear Systems

高斯消元法

回忆在线性代数中我们的做法，使用高斯消元，把矩阵化成上三角的形式，然后一个一个解出解并回代。

我们可以用一个迭代的思路来理解：我们先把最左列除了第一行的元素清零，然后，根据高斯消元法，我们就不再会对第一行、第一列操作了。这时，我们可以忽略它们，这就相当于是再去解一个新的、阶数减一的线性方程组了。迭代进行到阶数为 $1$ 的时候，解就出来了。

消元阶段的乘除法的数量为

$$\begin{align*} &\sum_{k=1}^{n-1}(n-k)(n-k-2)\\ &=\frac{n^3}{3}+\frac{n^2}{2}-\frac{5n}{6} \end{align*}$$

进行第 $k$ 次迭代时，需要对 $n-k$ 行做消去操作，每一行需要操作的元素有 $n-k+2$ 个。进行 $n-1$ 次操作后，矩阵阶数变为 $1$。

回代阶段的乘除法数量为

$$1+\sum_{i=1}^{n-1}(i+1)=\frac{n^2}{2}+\frac{n}{2}$$

考虑第 $i$ 次回代需要 $1$ 次除法和 $i$ 次乘法即可。

作业-3-2

P357-T8

Algorithm 6.1

大致步骤一致，只是在 Step4 部分有所区别：我们需要对 $j=1, \ldots n, j\ne i$ 操作。然后，我们在 Step8 阶段就可以返回所有的解：

P358-T11

根据上一题，我们可以得出，加减法的次数为

乘除法的次数为

为了便于计算，我们可以都只对乘除法的次数进行比较。高斯消元法需要的乘除法总数

故有

主元法 Pivoting Strategies

Partial pivoting/maximal column pivoting

第 $i$ 次迭代时， 选择第 $i$ 列中最小的 $p\gt i$，使得 $|a_{pi}|$ 是最大的。然后把第 $p$ 行和第 $i$ 行交换。这样可以保证除的时候除以一个较大的数。

然而，这个方法有一个问题，即除数大无法保证伸缩系数小。

时间复杂度： $O(n^2)$

Scaled Partial Pivoting/scaled-column pivoting

不仅仅看 $|a_{pi}|$，也考虑本行中的表现。

我们定义一个 scale factor $s_p=\max|a_{pj}|$。这表现了这一行中的伸缩系数分子中的最大值。

然后，我们寻找第 $i$ 列中最小的 $p\gt i$，使得 $|\frac{a_{pi}}{s_p}|$ 是最大(即最大的伸缩系数最小)的，然后交换。

但是，这么做也有个问题：每一次高斯消元的递归过程中，后面的每一行都发生了变化。如果我们每一行都这么做，会导致速度很慢。因此，注意 partial 这个词，我们的伸缩系数只在一开始计算，后面都不重新计算了。

时间复杂度： $O(n^2)$

Complete Pivoting/maximal pivoting

一个想法是，如果在 Scaled Partial Pivoting 中每次递归过程都重新计算 $s_p$，那么时间复杂度会到达立方。既然如此，不如也对列做交换。

对于 $n$ 阶的线性方程组，在第 $i$ 次递归时，会有一个 $n+1-i$ 阶子阵。在这个子阵中，找到最大的元素 $a_{pq}$，然后交换 $i$ 行与 $p$ 行，再交换 $i$ 列与 $q$ 列。

时间复杂度： $O(\frac{n^3}{3})$

矩阵分解 Matrix Factorization

让我们考虑每一次消元的等价操作：

令 $m_{i1}=\frac{a_{i1}}{a_{11}}(a_{11}\ne 0)$，则有第一次消元等价于将原系数矩阵左乘矩阵 $L_1$:

依次类推，第 $i$ 次消元等价于左乘 $L_i$:

既然如此，我们可以把每一个消元操作对应的矩阵先相乘。这样，可以得到一个下三角矩阵。

单位下三角矩阵乘与取逆都有封闭性。

这个下三角矩阵只与系数矩阵有关。因此，当系数矩阵一定时，这个方法效率很高。

现在，让我们重新整理一下符号。忘掉上面的 $L$，改取

则有

可以看到，我们把 $A$ 分解成了一个下三角矩阵 $L$ 和一个上三角矩阵 $U$ 的乘积。

如果我们在度量上做限制，让 $L$ 为一个单位下三角矩阵，那么这样的分解就是唯一的。

特殊矩阵 Special Types of Matrices

• 严格对角占优矩阵 Strictly Diagonally Dominant Matrix：每一行中，对角线元素绝对值大于本行其它元素的绝对值之和。严格对角占优矩阵是非奇异 nonsingular 的。且执行高斯消元法时不需要置换行或列，其解在舍入误差上是稳定的。

直觉上，严格对角占优矩阵类似于单位阵。因此，它会有比较好的性质。

• 正定矩阵 Positive Definite Matrix：对 $\forall x\ne \vec{0}, x^{-1}Ax>0$ 且 $A$ 是对称 symmetric 的.

  • 如果一个矩阵是正定的, 那么它的逆也是正定的.
  • 正定矩阵的对角线元素严格大于 $0$.

当给了矩阵的特殊性质，我们就可以针对性质给出特定的优化算法。

正定矩阵的优化算法

当 $A$ 是正定矩阵时，显然，它可以分解成 $B\Lambda B^T$ 的形式，其中 $B$ 为单位下三角矩阵， $\Lambda$ 为对角线元素均大于零的对角阵。另外，我们还有

其中 $\Lambda$ 是以 $u_{ii}$ 为对角元素的对角阵，而 $\tilde{U}$ 则因此变成了一个单位上三角矩阵。注意到 $L$ 为一个单位下三角矩阵，则只能有 $L=\tilde{U}^{T}$

又因为 $\Lambda$ 对角线元素均大于零，因此它可以被分解为 $\Lambda^{\frac{1}{2}}\Lambda^{\frac{1}{2}}$。那么，我们进行再整理，则有

使用 Choleski’s Method，可以快速地计算出 $\tilde{L}$，从而完成线性方程组的求解。

三对角矩阵的优化算法

Thomas 算法：

需要注意的是，一旦 $\exists\alpha_i =0$，则 Thomas 算法失效(但不意味着这个方程不可解)。

当 $A$ 是三对角矩阵，且它是 diagonally dominant 的(注意没有要求严格)，且

$$|b_1| >|c_1| >0, |b_n| > |a_n| >0, a_i\ne 0, c_i \ne 0$$

则 $A$ 非奇异，此时方程可用 Thomas 算法求解。

当 $A$ 是三对角矩阵且严格对角占优时，方程一定可用 Thomas 算法求解。

上述两种情况下，Thomas 算法是稳定的，因为所有的中间值会受主对角线元素的约束。

Thomas 算法的时间复杂度为 $O(n)$。

该算法的具体实现见作业题。

Homework-1

P397-T7

a)

我们采用一个递归的思路来解释 $LU$ 分解：

1. 在已知 $L$ 的前 $i-1$ 列， $U$ 的前 $i-1$ 行的情况下，求出 $U$ 的第 $i$ 行
2. 在已知 $L$ 的前 $i-1$ 列， $U$ 的前 $i$ 行的情况下，求出 $L$ 的第 $i$ 列

只要可以解决这两个递归，就能够完成 $LU$ 分解。

对于递归(1)，有(注意 $l_{ii}=1$ )：

这里引入了 $(i-1)(n+1-i)$ 次乘除法与 $(i-1)(n+1-i)$ 次加减法。

对于递归(2)，有：

这里引入了 $i(n-i)$ 次乘除法与 $(i-1)(n-i)$ 次加减法。

从 $i=1$ 一直递归到 $i=n$，就可以完成 $LU$ 分解了。

乘除法数为

加减法数为

b)

考虑 $\vec{y}$ 的前 $i-1$ 位已被求解，欲求解第 $i$ 位，有

这里引入了 $i-1$ 次乘除法和 $i-1$ 次加减法。

于是乘除法与加减法数都为

c)

与 (b) 中相反， $U\vec{x}=\vec{y}$ 的求解是从下往上的，且因为 $u_{ii}$ 可能不为 $1$，导致每一次都需要多进行一次除法，因此会每次引入 $i$ 次乘除法和 $i-1$ 次加减法

于是乘法数为

加减法数为

通过求和，可以解得总共需要的乘除法数量为

加减法数量为

这与高斯消元法需要的次数相同。

d)

乘除法总数为

加减法总数为

这个小题主要体现多次求解同个系数矩阵的线性方程组时 $LU$ 分解法的高效性。

Thomas 算法具体实现：

套用上面的 $LU$ 分解方法(虽然不是 ppt 中的 Crout 分解)，可以发现 $LU$ 分解的递归公式变得相当简单：

$$\begin{align*} u_{ii}&=a_{ii}-l_{ii-1}a_{i-1i}\\ u_{ii+1}&=a_{ii+1}\\ l_{ii}&=1\\ l_{i+1i}&=\frac{a_{i+1i}}{u_{ii}} \end{align*}$$

之后，我们再分别求解 $\vec{y}, \vec{x}$ 的递推式，有：

$$\begin{align*} y_{i+1}&=b_{i+1}-\frac{a_{i+1i}}{u_{ii}}y_i\\ x_{i-1}&=\frac{y_{i-1}-a_{i-1i}x_i}{u_{i-1i-1}} \end{align*}$$

注意到 $l$ 和 $u$ 的互相转换非常简单。因此，如果我们激进一些，完全可以只维护一个一维数组。我们记 $t_i=u_{ii}$，则有：

$$\begin{align*} t_{i+1}&=a_{i+1i+1}-\frac{a_{i+1i}a_{ii+1}}{t_i}\\ y_{i+1}&=b_{i+1}-\frac{a_{i+1i}}{t_i}y_i\\ x_{i-1}&=\frac{y_{i-1}-a_{i-1i}x_i}{t_{i-1}} \end{align*}$$

Read the proofs on P401-402

(i) 考虑 $A\vec{x}=0$，则有 $\vec{x}^{T}A\vec{x}=0$。根据正定性，有 $\vec{x}=\vec{0}$。因此，$A\vec{x}=0$ 只存在零解，故矩阵满秩，行列式不为零，是可逆的。

(ii) 这由正定性是显然的，我们只要取 $\vec{x_i}$ 的第 $i$ 位为 $1$，其它位为零，然后用 $\vec{x}^TA\vec{x}>0$ 的条件代入就可知， $A$ 的对角线元素必然严格大于零。

(iii) 我们同样采用构造的方法。我们任意取 $i, j$，使得 $x_i=1, x_j=-1$，且 $\vec{x}$ 的其它位为 $0$。那么，有(采用爱因斯坦求和约定)：

不妨设 $A_{ii}\ge A_{jj}$，则有

同理，我们把 $x_j$ 改成 $1$，就能得到 $A_{ii}>-A_{ij}$。

综上，对于任意的 $j$，均存在 $i\ne j$ 使得 $A_{ii}>|A_{ij}|$，证毕。

(iv) 令 $x_i=t, x_j=1, i\ne j$，其余位为 $0$，则有

左式作为一个一元二次方程，它的 $\Delta$ 应当小于零，即

化简即为所求。

P412-T17

(a) 只要计算行列式为 $0$ 即可。

(b) $\forall|\alpha| >1, |\beta| < 1$ 即可

(c) $\forall\alpha; \beta=1$

(d) 只需所有主子矩阵的行列式均为正。于是有：

得 $\beta < \frac{3}{2}\alpha \land \alpha>0$。

矩阵代数中的迭代方法 Iterative Techniques in Matrix Algebra

我们把 $A\vec{x}=\vec{b}$ 转换成迭代形式 $\vec{x}=T\vec{x}+\vec{c}$。为了能够迭代，我们先需要进行数学上的定义。

向量的范数

我们定义向量的范数 norm。范数需要满足:

• $k$ 阶范数的定义：

$$\|\vec{x}\|_p=\left(\sum_{i=1}^n|x_i|^p\right)^{\frac{1}{p}}$$

• 无穷阶范数的定义：

$$\|\vec{x}\|_{\infty}=\max_{1\le i\le n}|x_i|$$

欧几里得范数就是二阶范数。

对于某一个范数定义, 如果一个迭代序列 $\vec{x_k}$ 满足 $\|\vec{x_k}-\vec{x}\|$ 收敛于 $0$, 则说 $\vec{x_k}$ 收敛到 $\vec{x}$.

在有限维实向量空间中，不同范数之间总是“等价”的。也就是说，对于任意两个范数 $\|\cdot\|_a$ 和 $\|\cdot\|_b$，存在正的常数 $c$ 和 $C$，使得对所有 $x$ 都有

这说明：

• 如果一个迭代序列在范数 $\|\cdot\|_a$ 下收敛，那么在任意其他范数 $\|\cdot\|_b$ 下也收敛；
• 收敛的快慢（误差量级）在不同范数下最多相差一个常数因子，不影响“收敛阶”的讨论。

因此，在数值算法和迭代方法中，我们通常只选用一种范数来分析误差，而不必担心选哪一种范数会改变收敛性质。

矩阵的范数

一般而言，前面提到的性质对于范数已经足够了。然而，为了便于分析，我们往往要求矩阵的范数还满足如下条件：

如果没有一致性，我们对简单的 $AB$ 的误差分析都会变得很复杂，因为没办法从 $\|A\Delta B\|$ 到 $\|A\|\|\Delta B\|$ 了。

下面也介绍常见的范数：

• Fronenius 范数：

• $p$ 阶自然范数：

这个依赖于向量的定义似乎有些奇怪，但其实它表现了空间各个方向上矩阵算子最强的拉伸能力。

常用的自然范数有：

无穷自然范数的证明：

特征值与特征向量 Eigenvalues and Eigenvectors

谱半径 Spectral Radius

对于矩阵 $A$ ，其特征值可能是复数。因此，我们利用复数的幅值来刻画它们，其中最大的幅值也就是谱半径 $\rho(A)$ 。

对于一个 $n$ 阶矩阵，对于任意的自然范数，有 $\rho (A)\le \left\|A\right\|$

Prove:

$$\forall\text{ eigenvalue }\lambda, |\lambda|\cdot \left\|\mathbf{x}\right\|=\left\|\lambda \mathbf{x}\right\|=\left\|A \mathbf{x}\right\|\le \left\|A\right\|\cdot \left\|\mathbf{x}\right\|$$

如果一个 $n$ 阶矩阵满足对于任何 $i, j=1, 2, \ldots, n$，有 $\lim_{k\to\infty}(A^k)_{ij}=0$，则说这个矩阵是收敛的。

雅可比迭代法 Jacobi Iterative Method

雅可比迭代法的思路是构造一个便于求解的迭代式(系数矩阵只是个对角阵)，以提高迭代效率。

我们把一个矩阵 $A$ 分成对角阵 $D$ 与上三角阵 $-U$ 和下三角阵 $-L$ 的和，则有

这就出现了一个可迭代式(假设 $D$ 可逆)

我们定义雅可比迭代法的系数 $T_j=D^{-1}(L+U), \mathbf{c_j}=D^{-1}b$，则有迭代式

具体算法：

首先，我们合适地交换行或列，使得对于 $i=1, \ldots, n$， $a_{ii}\ne 0$。

然后，根据迭代式，每一次迭代，我们从 $i=1$ 直到 $n$，计算

$$x_i^{(k+1)}=\frac{b_i-\sum_{j\ne i \& j=1}^{n}a_{ij}x_j^{(k)}}{a_{ii}}$$

在这里，我们需要使用两个向量来维护 $\mathbf{x}^{(k)}, \mathbf{x}^{(k+1)}$。当不使用并行优化时，这么做显得有些浪费。一个想法是，我们始终只维护一个向量 $\mathbf{x}$ ，此时，在计算 $x_i$ 时，是使用了新的 $1\ldots i-1$ 个分量和其它旧的分量。这也就引出了高斯-赛德尔迭代法。

高斯-赛德尔迭代法 Gauss-Seidel Iterative Method

在直觉上，对于将 $A$ 分裂，然后一部分放到方程右边构造的迭代法，扔到方程右边的部分越多，收敛得就越快。

这通过极限思想容易想通，把 $A$ 拆成 $A$ 和 $0$ 后，一次迭代就收敛了。

当然，我们还需要保证迭代式易于求解的问题，因此一些拥有优良性质的矩阵就适合被留下来，把剩余部分扔到迭代右式里。

Gauss-Seidel Iterative Method 的思路就是解一个系数矩阵为下三角矩阵的迭代式。我们做如下拆分：

类似地，我们定义 $T_g=(D-L)^{-1}U, \mathbf{c_g}=(D-L)^{-1}\mathbf{b}$，则有迭代式

在具体实现上，很巧的是，这恰恰是我们在尝试对雅可比迭代法的空间占用上做优化的结果。

对于不同的问题，Jacobi 和 Gauss-Seidel 迭代法有不同的表现。有时前者好用后者不好用，反之的情况也会存在。但至少让我们关注一下收敛性的问题吧。

迭代方法的收敛性

首先介绍如下定理：

对于矩阵 $A$，下面五个命题等价：

1. $A$ 是收敛的。
2. 对于某个自然范数，有 $\lim_{n\to\infty}\left\|A\right\|=0$
3. 对于任意自然范数，有 $\lim_{n\to\infty}\left\|A\right\|=0$
4. $\rho (A)<1$
5. 对于任意的 $\mathbf{x}$，有 $\lim_{n\to\infty}A^n \mathbf{x}=\mathbf{0}$

拥有如上定理之后，我们可以得出一个结论：

对于任意的迭代式

$$\mathbf{x}^{(k)}=T\mathbf{x}^{(k-1)}+\mathbf{c}$$

只要满足 $\rho (T)<1$，迭代式就必然收敛到解

$$\mathbf{x}=T \mathbf{x}+\mathbf{c}$$

Proof:

我们将迭代式展开，则有

$$\mathbf{x}^{(k)}=T^k \mathbf{x}^{(0)}+(T^{k-1}+\ldots+T+I)\mathbf{c}$$

由于 $\rho(T)<1$，有对于任意的 $\mathbf{x}$， $\lim_{n\to\infty}T^n \mathbf{x}=\mathbf{0}$。因此，当我们取 $k$ 无限大时，就有 $T^k \mathbf{x}^0\to \mathbf{0}$。

于是，我们有

$$\lim_{k\to\infty}\mathbf{x}^{(k)}=(I-T)^{-1}\mathbf{c}$$

而这恰好就是 $\mathbf{x}=T \mathbf{x}+\mathbf{c}$ 的解。

回忆前文，我们知道，对于任意的自然范数 $\left\|\cdot\right\|$，有 $\rho (T)\le \left\|T\right\|$。这说明， $\left\|T\right\|<1$ 是比 $\rho (T)<1$ 更强的条件。当满足这个更强的条件时，有如下结论：

迭代式

$$\mathbf{x}^{(k)}=T\mathbf{x}^{(k-1)}+\mathbf{c}\quad(\text{given }\forall \text{ natural norm}\left\|T\right\|<1)$$

收敛，且其中间解与收敛终点的误差控制在：

$$\begin{align} \left\|\mathbf{x}-\mathbf{x}^{(k)}\right\|&\le \left\|T\right\|^k \left\|\mathbf{x}-\mathbf{x}^{(0)}\right\|\\ \left\|\mathbf{x} -\mathbf{x}^{(k)}\right\|&\le \frac{\left\|T\right\|^k}{1-\left\|T\right\|}\left\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\right\| \end{align}$$

对于一个严格对角占优的矩阵 $A$，不管迭代起点如何选取，Jacobi 和 Gauss-Seidel 法都会得到收敛到唯一解的迭代序列。

Relaxation Methods

回忆第一个 research topic 中我最终对不动点迭代的方法，可以发现一个特殊的视角：

我们可以把不动点迭代写成

$$\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+ \mathbf{\Delta }^{(k+1)}$$

其中 $\mathbf{\Delta }^{(k+1)}$ 是一个有关于 $\mathbf{x}^{(k)}$ 的函数，且当 $\mathbf{x}^{(k)}$ 为不动点时取值为 $\mathbf{0}$。这个函数事实上表现了每一步迭代的步长。

接着这个视角，我们可以把 Gauss-Seidel 迭代法写成如下的形式：

特殊地，我们记

并称之为残差 residual。

这样，就有

回忆一下我们上面提到的视角，我们很容易想到一件事，这个迭代的步长也许可以加一些系数来调整。这也就是所谓的 relaxation method：

根据 $\omega$ 的范围不同，我们把 relaxation method 分为

• $0<\omega <1$ Under-Relaxation Methods
• $\omega =1$ Gauss-Seidel Methods
• $\omega >1$ Successive Over-Relaxation Methods