研一课程笔记-数值分析1-3

绪论与预备知识

概论

衡量数值算法的几个标准

衡量算法的优劣有两个标准，一是要有可靠的理论基础，包括收敛性、数值稳定性和算法精度等；二是要有较好的计算复杂度。

• 数值稳定性

  • 数值稳定性即算法对舍入误差的敏感性。舍入误差对计算结果的精确性影响小的算法，具有较好的数值稳定性；反之，算法的数值稳定性差。

• 算法精度

  • 数值算法的精度衡量了近似解对精确解的近似程度。算法的精度越高并不一定导致近似解的误差越小。精度的具体定义依赖于所考虑的问题类型。在偏微分方程的数值算法中我们一般考虑收敛精度，它指的是算法的收敛速率。

• 计算复杂度

  • 计算复杂度指的是通过数值计算解决问题的困难程度。最常见的是时间复杂度和空间复杂度

• 收敛性

  • 算法的收敛性一般指数值解逼近精确解的性质，而对于不同的问题，其定义也是有区别的。例如在迭代算法中，收敛性指的是近似解能够随着迭代过程趋向于精确解，它通常与迭代矩阵的谱半径相关。在偏微分方程数值算法中，收敛性指的是随网格尺寸等参数变化的数值解的收敛性质，它与算法的收敛精度密切相关。

误差理论

误差按照它们的来源可分为以下四种

1. 模型误差：反映实际问题有关量之间关系的计算公式，即数学模型，通常只是近似的。由此产生的数学模型的解与实际问题的解之间的误差称为模型误差。

2. 观测误差：数学模型中包含的某些参数（如时间、长度、电压等等）往往通过观测而获得。由观测得到的数据与实际的数据之间是有误差的。这种误差称为观测误差。

3. 截断误差求解数学模型所用的数值计算方法如果是一种近似的方法，那么只能得到数学模型的近似解，由此产生的误差称为截继误差或方法误差。例如，由Taylor (泰勒)公式，函数 $f(x)$ 可表示为

$$f(x)=f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n+\frac{f^{(n+1)}(\theta x)}{(n+1)!}x^{n+1},(0<\theta<1)$$

近似即：（此近似公式的误差就是截断误差）

$$f(x)\approx f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n$$

4. 舍入误差：由于计算机的字长有限，参加运算的数据以及运算结果在计算机上存放会产生误差

数值计算的误差度量

对于准确值$x$、近似值 $a$ 和误差限 $\epsilon$ ：

• $e=x-a$ ，称 $e$ 为近似值 $a$ 的绝对误差，简称误差。

• 如果 $|e|$ 的一个上界已知，记为 $\epsilon$ ，即 $|e|\leq \epsilon$ 则称 $\epsilon$ 为近似值 $a$ 的绝对误差限或绝对误差界，简称误差限或误差界。

• $e_r=\frac ex=\frac{x-a}x$ ，称 $e_r$ 为近似值 $a$ 的相对误差。由于 $x$ 未知，实际上总把 $\frac e a$作为 $a$ 的相对误差，并且也记为$e_r=\frac ea=\frac{x-a}a$。一般用百分比表示

• $|e_r|$ 的上界，即 $\epsilon_r=\frac{\epsilon}{|a|}$ 称为近似值 $a$ 的相对误差限或相对误差界。

• 有效数字：由准确值经过四舍五入得到的近似值（经过对于某一绝对误差限舍入后），从它的末位数字到第一位非零数字都是有效数字

  • 严格来说：如果 $a$ 的绝对误差限是它的某一位的半个单位，并且从该位到它的第一位非零数字共有 $n$ 位，则称用 $a$ 近似 $x$ 时具有 $n$ 位有效数字。

PS. $a-\epsilon\leq x\leq a+\epsilon$，即 $x=a\pm\epsilon.$

函数值的误差估计

给定多元函数$A=f\left(x_1,x_2,\cdots,x_n\right)$ ，且 $x_1^*,x_2^*,\cdots,x_n^*$ 为 $x_1,x_2,\cdots,x_n$ 近似值，于是可求 $A$ 的近似值 $A^*=f\left(x_1^*,\cdots,x_n^*\right)$

绝对误差

在点 $x=\left(x_1,x_2,\cdots,x_n\right)$ 进行泰勒展开，由

$$f(x_{1}, x_{2}, \ldots, x_{n}) = f(a_{1}, a_{2}, \ldots, a_{n}) + \sum_{i=1}^{n} \frac{\partial f}{\partial x_{i}} (x_{i} - a_{i}) + \frac{1}{2!} \sum_{i,j=1}^{n} \frac{\partial^{2} f}{\partial x_{i} \partial x_{j}} (x_{i} - a_{i})(x_{j} - a_{j}) + \cdots$$

其中， $\frac{\partial f}{\partial x_{i}}$ 表示函数 $f$ 对变量 $x_{i}$ 的偏导数，$\frac{\partial^{2} f}{\partial x_{i} \partial x_{j}}$ 表示函数$f$ 对变量 $x_{i}$ 和 $x_{j}$ 的混合偏导数，依此类推。所以：

$$\begin{aligned}A^{*}-A&=f\left(x_{1}^{*},\cdots,x_{n}^{*}\right)-f\left(x_{1},x_{2},\cdots,x_{n}\right)\\&\approx\sum_{j=1}^n\frac{\partial f(x)}{\partial x_j}\left(x_j^*-x_j\right)\end{aligned}$$

PS.这个等式也可f以用中值定理直接得出

由 $e(x_i)=x_j^*-x_j$ ，也即：

$$e(A)\approx\sum_{j=1}^n\frac{\partial f(x)}{\partial x_j}e\left(x_j\right),|e(A)|\leqslant\sum_{j=1}^n\left|\frac{\partial f(x)}{\partial x_j}\right||e\left(x_j\right)|$$

PS.

相对误差

$$e_r(A)=\frac{e(A)}{A}\approx\sum_{j=1}^n\frac{x_j}{f(x)}\frac{\partial f(x)}{\partial x_j}e_r\left(x_j\right)$$

• $e_r (x_j) = \frac{e(x_j)}{x_j}$ ，即 $e(x_j)= e_r (x_j)\cdot {x_j}$
• 因子 $\frac{x_j}{f(x)}\frac{\partial f(x)}{\partial x_j}$ 反映 $x_j$ 相对误差 $e_r (x_j)$ 对相对误差 $e_r(A)$ 影响的程度

函数值误差的代数运算

• $f(x,y)=xy$
  • $e(xy)\approx ye(x)+xe(y)$
  • $e_r(xy)\approx e_r(x)+e_r(y)$
• $f(x,y)=x/y$
  • $e(x/y)\approx\frac{1}{y}e(x)-\frac{x}{y^{2}}e(y)$
  • $e_{r}(x/y)\approx e_{r}(x)-e_{r}(y)$
• $f(x,y)=x\pm y$
  • $e(x\pm y)\approx e(x)\pm e(y)$
  • $e_{r}(x\pm y)\approx\frac{x}{x\pm y}e_{r}(x)\pm\frac{y}{x\pm y}e_{r}(y)$ ， $x\pm y\neq0$
• $f(x)=\sqrt{x}$
  • $e_r(\sqrt{x})\approx\frac{1}{2}e_r(x)$
• $f(x)=x^{n}$
  • $e_{r} (x^{n})\approx ne_{r}(x)$

几点注意：

• 对于 $e_{r}(x- y)\approx\frac{x}{x- y}e_{r}(x)-\frac{y}{x- y}e_{r}(y)$ ，如果 $x$ 与 $y$ 非常接近，会出现很大计算误差
• 同样的，对于 $e(x/y)\approx\frac{1}{y}e(x)-\frac{x}{y^{2}}e(y)$ ，若 $y\to 0$，同理

算法在数值计算中应注意的几个问题

• 要有数值稳定性。即能控制舍入误差的传播
• 两数相加要防止较小的数加不到较大的数中所引起的严重后果
• 避免两个相近的近似值相减，以免严重损失有效数字
• 除法运算中，要尽量避免除数的绝对值远小于被除数的绝对值

范数

向量范数

定义向量大小的量，又称为向量的模。是定义在 $R^n$ 上的实值函数 $\|\cdot\|$，满足：

• $\|x\|\geq 0$ 且仅 $x=0,\|x\|=0$
• $\|kx\|=|k|\cdot\|x\|,k\in R$
• $\|x+y\|\leq\|x\|+\|y\|$

证明是向量范数，只需要满足上面三个条件

三种向量范数($p-$范数)

对 $R^n$ 中的任一向量 $x=(x_{1},x_{2},\cdots,x_{n})^{\mathrm{T}}$

$$\|x\|_{1}=\sum_{i=1}^{n}|x_{i}|\\\|x\|_{2}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}\\\|x\|_{\infty}=\max_{1\leqslant i\leqslant n}|x_{j}|$$

范数等价性定理

设 $\|\cdot\|_\alpha,\|\cdot\|_\beta$ 是 $R^n$ 上的任意两种向量范数，则存在与向量 $x$ 无关的常数 $m$ 和 $M(0 < m < M)$， 使下列关系成立

$$m\|x\|_\alpha\leqslant\|x\|_\beta\leqslant M\|x\|_\alpha,\quad\forall x\in\mathbb{R}^n$$

此外：

摘自 张绍飞，赵迪. 矩阵论教程[M]. 机械工业出版社，2012. p98

矩阵范数

用于定义矩阵“大小”的量，又称为矩阵的模，定义在 $R^{n\times n}$ 上的实值函数 $\|\cdot\|$，满足对任意矩阵 $A,B$：

• $\|A\|\geq 0$ 且仅 $A=O,\|A\|=0$
• $\|kA\|=|k|\cdot\|A\|,k\in \mathbb R$
• $\|A+B\|\leq\|A\|+\|B\|$
• $\|AB\|\leq\|A\|\cdot\|B\|$

矩阵范数相容性

对于给定的向量范数 $\|\cdot\|$ 和矩阵范数 $\|\cdot\|$ ， 如果对任一个 $x\in \mathbb R^n$ 和任一个 $A\in \mathbb R^{n\times n}$， 满足 $\|Ax\|\leqslant\|A\|\|x\|$ 则称所给的矩阵范数与向量范数是相容的

当在同一个问题中需要同时使用矩阵范数和向量范数时，这两种范数应当是相容的

矩阵的算子范数定理

设在 $R^n$ 中给定了一种向量范数，对任一矩阵 $A\in\mathbb{R}^{n\times n}$， 令：

$$\|A\|=\max_{\|x\|=1}\|Ax\| \tag{1.3.1}$$

则由式1.3.1定义的 $\|\cdot\|$是一种矩阵范数，并且它与所给定的向量范数相容

称式1.3.1所定义的矩阵范数为从属于所给定向量范数的矩阵范数，又称为矩阵的算子范数

$\max_{\|x\|=1}\|Ax\|$ 表示矩阵 $A$ 作用在单位范数向量 $x$ 上所得到的结果的最大范数。具体来说：

• 对于给定的矩阵 $A$，我们考虑所有单位范数向量 $x$，即满足 $\|x\|=1$ 的向量。
• 我们将矩阵 $A$ 作用在每一个单位范数向量 $x$ 上，得到向量 $Ax$。
• 然后计算每个 $Ax$ 的范数 $\|Ax\|$。
• 最终，我们找到使得 $\|Ax\|$ 最大的单位范数向量 $x$，即求解 $\max_{\|x\|=1}\|Ax\|$。

证明为什么是一种矩阵范数：

设给定的向量范数为 $\|\cdot\|_p$，则从属于向量范数 $\|\cdot\|_p$ 的矩阵范数仍记为 $\|\cdot\|_p$，即

$$\|A\|_p=\max_{\|x\|_p=1}\|Ax\|_p$$

矩阵 $p-$范数

设 $A=[a_{ij}]\in\mathbb{R}^{n\times n}$，则

$$\begin{gathered} \|A\|_{1}=\max_{1\leqslant j\leqslant n}\sum_{i=1}^{n}|a_{ij}| \\ \|A\|_{2}=\sqrt{\lambda_{\max}\left(A^{\mathrm{T}}A\right)} \\ \|A\|_{\infty}=\max_{1\leqslant i\leqslant n}\sum_{j=1}^{n}|a_{ij}| \end{gathered}$$

其中 $\lambda_{\mathrm{max}} (A^{\mathrm{T}}A)$表示矩阵 $A^{\mathrm{T}}A$ 的最大特征值。

矩阵范数 $\|A\|_{1},\|A\|_{2},\|A\|_{\infty}$ 又分别称为矩阵的列范数、谱范数和行范数

佛罗贝尼乌斯范数

另外还有一种常用的矩阵范数，就是Frobenius(佛罗贝尼乌斯) 范数，又称为Euclid 范数，

$$\|A\|_\mathrm{F}=\sqrt{\sum_{i,j=1}^na_{ij}^2}$$

其中 $A=[a_{ij}]\in\mathbb{R}^{n\times n}$，又可以记作 $\|A\|_\mathrm{E}$

可以证明,$\|\cdot\|_\mathrm{F}$ 与向量范数 $\|\cdot\|_{2}$ 相容，即

$$\|Ax\|_2\leqslant\|A\|_\text{F}\|x\|_2,\quad A\in\mathbb{R}^{n\times n},x\in\mathbb{R}^n$$

但 $\|\cdot \|_\mathrm{F}$ 不从属于任何向量范数，即不是算子范数

算子范数的估计

单位矩阵 $I$ 的任何一种算子范数都有

$$\|I\|=\max\limits_{\|x\|=1}\|Ix\|=1$$

---

设矩阵 $A \in R^{n\times n}$ 的某种范数 $\left\|A\right\|<1$ ，则 $I\pm A$ 为非奇异矩阵，并且当该种范数为算子范数时，下式成立

$$\left\|(I\pm A)^{-1}\right\|\leqslant\frac1{1-\|A\|}$$

若 $I\pm A$ 奇异，即 $(I\pm A)\cdot x=0$，展开并移项，即 $x=\mp Ax$，取范数，即 $\|x\|=\|Ax\| \leq \|A\|\cdot \|x\|$

即 $\|A\|\geq1$，矛盾

当然这里还涉及到一个点，即矩阵范数相容性，即满足 $\|Ax\|\leqslant\|A\|\|x\|$

此外也可以通过这种方式来证明：

取

$$e_i=\left(\begin{array}{c}0\\1\\0\end{array}\right)$$

则：

$$x\cdot e_i=\begin{bmatrix} 0 & x_1 & 0 \\ 0 & x_2 & 0 \\ 0 & x_3 & 0 \end{bmatrix}$$

取范数，即 $\|x\|=\|x\cdot e_i\|$，相当于我们构建了一个向量范数与矩阵范数的等式

那么 $\|Ax\|=\|Ax\cdot e_i\|\leq\|A\|\cdot\|x\cdot e_i\|=\|A\|\cdot\|x\|$

由 $(I-A)(I-A)^{-1}=I$ ，则 $(I-A)^{-1}-A\cdot (I-A)^{-1}=I$，即 $(I-A)^{-1}=A\cdot (I-A)^{-1}+I$

两边同取算子范数：$\|(I-A)^{-1}\|=\|A\cdot (I-A)^{-1}\|+\|I\|\leq\|A\|\cdot\|(I-A)^{-1}\|+\|I\|$

移项，即 $(1-\|A\|)\|(I-A)^{-1}\|\leq\|I\|=1$，而 $\|A\|<1$，故上式得证

PS.一个方阵是奇异的，如果它的行列式为零。非奇异矩阵则是行列式非零的方阵。

1. 逆矩阵：奇异矩阵是不可逆的，因为逆矩阵的存在要求矩阵的行列式不为零。因此，奇异矩阵没有逆矩阵。

2. 秩：奇异矩阵的秩小于其阶数。一个 $n \times n$ 的奇异矩阵的秩最多为 $n-1$。

3. 零空间：奇异矩阵的零空间（也称为核）不仅包含零向量，还包含其他非零向量。这是因为奇异矩阵将某些向量映射到零向量。

4. 特征值：奇异矩阵的特征值中至少有一个为零。这是因为特征值是行列式与矩阵迹的根，而奇异矩阵的行列式为零。

5. 解的存在性：对于线性方程组 $Ax = b$，其中 $A$ 是奇异矩阵，解可能不存在，或者有无穷多个解。这取决于向量 $b$ 是否位于矩阵 $A$ 的列空间中。

线性方程组的解法

主要讨论如下形式的线性方程组的求解问题.

$$\left.\left\{\begin{array}{c}a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1,\\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2,\\\vdots\\a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n,\end{array}\right.\right.\tag{2.1.1}$$

也即：

$$Ax=b$$

其中

$$\left.A=\left(\begin{array}{ccc}a_{11}&\cdots&a_{1n}\\\vdots&&\vdots\\a_{n1}&\cdots&a_{nn}\end{array}\right.\right)\in\mathbb{R}^{n\times n},\quad x=\left(\begin{array}{c}x_1\\\vdots\\x_n\end{array}\right)\in\mathbb{R}^n,\quad b=\left(\begin{array}{c}b_1\\\vdots\\b_n\end{array}\right)\in\mathbb{R}^n.$$

设系数矩阵 $A$ 非奇异，即 $det(A) \neq 0$，则方程组有唯一解向量 $x$.

Gauss 消去法

由消元和回代两部分组成

消元即是对增广矩阵 $[A,b]$ 做初等行变换，经过 $n-1$ 次消元得到

$$\left.\left(A^{(n)},b^{(n)}\right)=\left(\begin{array}{ccccc}a_{11}^{(1)}&\cdots&\cdots&a_{1n}^{(1)}&b_{1}^{(1)}\\&a_{22}^{(2)}&\cdots&a_{2n}^{(2)}&b_{2}^{(2)}\\&&\ddots&\vdots&\vdots\\O&&&a_{nn}^{(n)}&b_{n}^{(n)}\end{array}\right.\right)$$

即线性方程组：

$$\begin{aligned}\begin{cases}a_{11}^{(1)}x_1+a_{12}^{(1)}x_2+\cdots+a_{1n}^{(1)}x_n=b_1^{(1)},\\a_{22}^{(2)}x_2+\cdots+a_{2n}^{(2)}x_n=b_2^{(2)},\\\vdots\\a_{nn}^{(n)}x_n=b_n^{(n)}.\end{cases}\end{aligned}$$

通过逐步回代，依次求出 $x_{n-1},x_{n-2},\ldots,x_1$

顺序Gauss 消去法

消元过程的第 $k$ 步，对矩阵需做 $(n-k)^2$ 次乘法运算及 $(n - k)$ 次除法运算，对右端向量需作 $(n - k)$ 次乘法运算，所以消去过程总的乘除法运算工作量为

$$\sum_{k=1}^{n-1}(n-k)^2+\sum_{k=1}^{n-1}(n-k)+\sum_{k=1}^{n-1}(n-k)=\frac13n^3+\frac12n^2-\frac56n.$$

回代过程中，计算每个 $x_k$ 需作 $(n - k + 1)$ 次乘除法运算，其工作量为

$$\sum_{k=1}^n(n-k+1)=\frac12n(n+1).$$

算法运行条件

要使得顺序消去法可运行，必须使得前 $n - 1$ 个主元素 $a^{(k)}_{kk} (k = 1，2， \dots，n - 1)$ 均不为零，其充要条件是系数矩阵 $A$ 的前 $n-1$ 个顺序主子式不为零，即

$$\left.D_k=\left|\begin{array}{ccc}a_{11}^{(1)}&\cdots&a_{1k}^{(1)}\\\vdots&&\vdots\\a_{k1}^{(1)}&\cdots&a_{kk}^{(1)}\end{array}\right.\right|\neq0,\quad k=1,2,\cdots,n-1.$$

为什么运行条件是主元素不为零充要是系数主子式不为零，这个其实也可以用LU分解的思路来思考

当主元素 $a^{(k)}_{kk}$ 均不为零时，我们可以考虑使用LU分解的思路来证明对应系数矩阵的顺序主子式不为零。

假设我们有一个 $n \times n$ 的方阵 $A$，其顺序主子式都不为零，即 $|A_k| \neq 0$，其中 $A_k$ 表示 $A$ 的第 $k$ 阶顺序主子式。

我们知道LU分解可以写为 $A = LU$，其中 $L$ 是单位下三角矩阵，$U$ 是上三角矩阵。我们可以将 $A$ 分解为：

$$A = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{bmatrix}$$

我们可以看到，对于LU分解中的 $U$ 矩阵，其对角线上的元素 $u_{kk}$ 对应于原矩阵 $A$ 的主元素 $a^{(k)}_{kk}$。由于主元素均不为零，因此LU分解中的 $U$ 矩阵的对角线上的元素 $u_{kk}$ 也均不为零。

现在我们来看LU分解中的 $L$ 矩阵。由LU分解的性质可知，$L$ 的对角线上的元素均为1，而且 $L$ 的第 $k$ 行第 $j$ 列的元素 $l_{kj}$ 是原矩阵 $A$ 的第 $j$ 阶顺序主子式与 $U$ 的第 $k$ 阶顺序主子式的商。由于 $A$ 的顺序主子式都不为零，因此 $L$ 的元素 $l_{kj}$ 也不为零。

一个问题

前面我们提到，除法运算中，要尽量避免除数的绝对值远小于被除数的绝对值

如果对角线元素 $a_{ii}$ 很小（导致 $|m_{ik}|$ 很大），就会导致很大的误差

列主元素Gauss 消去法

第 $k$ 次消元之前，通过行变换将 $a^{(k)}_{ik} (i =k，k + 1,\dots ，n)$ 中绝对值最大的元素交换到第 $k$ 行的主对角线位置上，然后进行消元计算，即：

其他过程同顺序法，算法要求各个列的主元素不为零

对于列主元素交换，相当于多乘了一个初等行变换矩阵，所以最终还是可以化为一个LU分解的形式

比如乘以

$$P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

交换一二行

LU分解及其与Gauss消去的关系

考试不考，仅供拓展思路

对于 $Ax=b$ 将 $A$ 分解为 $A=LU$，其中 $L$ 是下三角矩阵， $U$ 是上三角矩阵，则可将原方程化为：

$$Ly=b,\,\,Ux=y$$

先由 $Ly=b$ 解出 $y$，再由 $Ux=y$ 解出 $x$

而Gauss消去也可以看作一个LU分解的过程，即：

$$L^{(n-1)}L^{(n-2)}\cdots L^1 \cdot A=A^{(n)}$$

其中 $L^1$ 对应于将第一列非主对角线元素都化为零的过程，即

$$L^1=\begin{bmatrix} 1 & 0 & 0 & \cdots &0\\ l_1 & 1 & 0 & \cdots &0\\ l_2 & 0 & 1 & \cdots &0\\ \vdots & \vdots & \vdots &&0\\ l_{n-1} & 0 & 0 & \cdots &1\\ \end{bmatrix}$$

其中 $l_i$ 对应 $-m_{i1}$ 。以此类推

不妨设 $L=(L^{(n-1)}L^{(n-2)}\cdots L^1)^{-1},U=A^{(n)}$，显然，$L$ 是下三角矩阵，$A$ 为上三角矩阵，即为 $A=LU$

矩阵条件数与病态方程组

设$A$、$\Delta A\in\mathbb{R}^n\times n$，$b$、$\Delta b\in\mathbb{R}^n$，$A$非奇异，$b\neq0$，$x$ 是方程组 $Ax=b$ 的解向量。若 $\|\Delta A\|<\frac1{\|{A}^{-1}\|}$ ，则有

• 方程组

  $$(A+\Delta A)(x+\Delta x)=b+\Delta b$$

  有唯一解 $x+\Delta x$

• 下列估计式成立：

  $$\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\frac{\parallel A\parallel\parallel A^{-1}\parallel}{1-\parallel A^{-1}\parallel\parallel\Delta A\parallel}(\frac{\parallel\Delta A\parallel}{\parallel A\parallel}+\frac{\parallel\Delta b\parallel}{\parallel b\parallel})$$

对于第一点，我们需要证明 $A+\Delta A$ 非奇异，证明关键其实是该定理：

设矩阵 $A \in R^{n\times n}$ 的某种范数 $\left\|A\right\|<1$ ，则 $I\pm A$ 为非奇异矩阵，并且当该种范数为算子范数时，下式成立

$$\left\|(I\pm A)^{-1}\right\|\leqslant\frac1{1-\|A\|}$$

由 $A+\Delta A=A(I+A^{-1}\Delta A)$，其中 $A$ 非奇异，即证 $I+A^{-1}\Delta A$ 非奇异，由上定理，即证 $\|A^{-1}\Delta A\|<1$

而 $\parallel A^{-1}\Delta A\parallel\leqslant\parallel A^{-1}\parallel\parallel\Delta A\parallel<1$，故得证

对于第二点，对原方程组进行展开，并代入 $Ax=b$

$$Ax+A\cdot \Delta x+\Delta A \cdot x+\Delta A\cdot \Delta x=b+\Delta b\\ A\cdot \Delta x+\Delta A \cdot x+\Delta A\cdot \Delta x=\Delta b\\ \Delta{x}={A}^{-1}\Delta{b}-{A}^{-1}\Delta{A}{x}-{A}^{-1}\Delta{A}\Delta{x}$$

两边同时取范数，由矩阵范数相容性，即 $\|b\|=\|Ax\|\leqslant\|A\|\|x\|$

$$\begin{aligned} \|\Delta x\| &= \|A^{-1}\Delta b\|+\|A^{-1}\Delta A\cdot x\|+\|A^{-1}\Delta{A}\Delta{x}\|\\ &\leq \parallel A^{-1}\parallel\parallel\Delta b\parallel+\parallel A^{-1}\parallel\parallel\Delta A\parallel\parallel x\parallel+\parallel A^{-1}\parallel\parallel\Delta A\parallel\parallel\Delta x\parallel \end{aligned}$$

即

$$(1-\parallel A^{-1}\parallel\parallel\Delta A\parallel)\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\parallel A\parallel\parallel A^{-1}\parallel\frac{\parallel\Delta b\parallel}{\parallel b\parallel}+\parallel A^{-1}\parallel\parallel\Delta A\parallel$$

由第一点 $\parallel A^{-1}\parallel\parallel\Delta A\parallel<1$，上式得证

一个细节就是这里用到的都是算子范数，使得 矩阵范数与向量范数相容，从而满足范数相容性

对于上面定理第二点的式子，我们可以改写成：

$$\frac{\parallel\Delta x\parallel}{\parallel x\parallel}\leqslant\frac{\parallel A\parallel\parallel A^{-1}\parallel}{1-\parallel A\parallel\parallel A^{-1}\parallel\frac{\parallel\Delta A\parallel}{\parallel A\parallel}}(\frac{\parallel\Delta A\parallel}{\parallel A\parallel}+\frac{\parallel\Delta b\parallel}{\parallel b\parallel})$$

那么 $\parallel A\parallel\parallel A^{-1}\|$ 越小，$A,b$ 的误差对解向量的相对误差影响就越小，反之则大。

条件数与病态

关于病态方程组求解

对于病态方程组的求解：

• 用更高精度的运算
• 平衡方法：$A$ 元素数量级差别很大时进行行平衡或列平衡
• 残差矫正法：病态但不特别严重时的一种缓解病态方法

迭代法

迭代法是目前求解大规模稀疏线性方程组的主要方法之一.

迭代法的一般形式及其收敛性

对于线性方程组，将系数矩阵分解为两个矩阵的差，即 $A=N-P$ ，代入 $Ax=b$，即：

$$Nx=Px+b\\x=N^{-1}Px+N^{-1}b$$

记 $G=N^{-1}P,\quad d=N^{-1}b$ 代入，得 $x=Gx+d$ ，取 $x^{(0)}\in\mathbb{R}^{n\times n}$ 作为初始近似解，使用递推式

$$x^{(k+1)}=Gx^{(k)}+d\quad(k=0,1,\cdots)$$

产生一个向量序列，当 $k$ 足够大时， $x^k$ 可以看作原方程的近似解

序列收敛性

只有在迭代法收敛的情况下，用它所产生的向量序列 $\{x^{(k)}\}$ 中的向量作为方程组的近似解才有意义，而且，k 越大，$x^{(k)}$ 作为方程组的解就越精确。关于向量序列收敛定义如下：

设有向量序列

$$x^{(k)}=\left({x_1}^{(k)},{x_2}^{(k)},\cdots,{x_n}^{(k)}\right)^\top\quad(k=0,1,\cdots)$$

如果存在常向量$x^*=(x_1^*,x_2^*,\cdots,x_n^*)^\top$，使得

$$\lim\limits_{k\to\infty}x_i^{(k)}=x_i^*\quad(i=1,2,\cdots,n)$$

则称向量序列$\{x^(k)\}$收敛于常向量$x^*$，记为

$$\lim\limits_{k\to\infty}x^{(k)}=x^*$$

相关定理

设有向量序列$\{x^{(k)}\}$和常向量 $x^*$，如果对某种范数，有

$$\lim\limits_{k\to\infty}\left\|x^{(k)}-x^*\right\|=0$$

则必有

$$\lim\limits_{k\to\infty}x^{(k)}=x^*$$

由向量范数等价性，原条件即：

$$\lim_{k\to\infty}\max_{1\leqslant i\leqslant n}\left|x_i^{(k)}-x_i^*\right|=0$$

一个后面用到的等式

$$x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)$$

要证明这个公式，我们从Jacobi迭代法的定义出发。Jacobi迭代法用于求解线性方程组 $Ax = b$，迭代公式为：

$$x^{(k+1)} = D^{-1}(b - (L + U)x^{(k)})$$

其中 $A = D + L + U$，$D$ 是 $A$ 的对角矩阵，$L$ 是 $A$ 的严格下三角矩阵，$U$ 是 $A$ 的严格上三角矩阵。

我们要证明：$x^{(k)} - x^* = \sum_{i=k}^\infty (x^{(i)} - x^{(i+1)})$

首先，我们知道对于任何迭代步 $k$，有以下关系：

$$x^{(k)} - x^* = (x^{(k)} - x^{(k+1)}) + (x^{(k+1)} - x^*)$$

利用这一点，我们可以递推：

$$x^{(k+1)} - x^* = (x^{(k+1)} - x^{(k+2)}) + (x^{(k+2)} - x^*)$$

这样继续下去，我们得到：

$$x^{(k)} - x^* = (x^{(k)} - x^{(k+1)}) + (x^{(k+1)} - x^{(k+2)}) + (x^{(k+2)} - x^{(k+3)}) + \cdots$$

在极限的意义下：

$$x^{(k)} - x^* = \sum_{i=k}^\infty (x^{(i)} - x^{(i+1)})$$

这个公式成立的关键在于迭代过程的收敛性，即 $x^{(i)} \to x^*$ 当 $i \to \infty$。因此，尾项 $(x^{(i)} - x^*)$ 在极限下趋于零。

谱半径与收敛性

设 $n\times n$ 矩阵 $G$ 的特征值是 $\lambda_1,\lambda_2,\cdots,\lambda_n$ ,称

$$\rho(G)=\max_{1\leqslant i\leqslant n}|\lambda_i|$$

为矩阵 $G$ 的谱半径.

谱半径与收敛性的定理

• 定理：任意一个矩阵的谱半径与其范数是有关系的，即任一矩阵$A$的谱半径不大于$A$的任一与某一向量范数相容的矩阵范数，即

$$\rho(A)\leqslant\|A\|.$$

按照谱半径的定义 $\rho(A)=\max_{1\leq i\leq n}|\lambda_i|$， 其中，$\lambda_i(i=1,2,\cdots,n)$为 $n$ 阶方阵 $A$ 的特征值

设 $x$ 为对应于 $\lambda$ 的 $A$ 的特征向量，则有 $\lambda x=Ax$，两边取范数可得

$$|\lambda|\|x\|\leqslant\|A\|\|x\|.$$

因为 $x$ 为非零向量，$\|x\|\neq0$，故有$|\lambda|\leqslant\|A\|.$

上述这个不等式对$A$的任何特征值均成立，即可得 $\rho(A)\leqslant\|A\|.$

• 对上定理，利用矩阵的Jordan 标准型得：设 $A$ 是任意 $n$ 阶矩阵, 则 $A$ 的 $k$ 次幂 $A^k\to 0$ ( 当 $k \to +\infty$) 的充要条件为谱半径 $\rho(A)<1$

Jordan标准型是一种特殊形式的矩阵，通常用于简化矩阵的分析和计算特征值、特征向量等问题。对于一个给定的方阵，其Jordan标准型可以通过相似变换得到。以下是关于Jordan标准型的一些重要概念：

1. 定义：对于一个 $n \times n$ 的方阵 $A$，存在一个可逆矩阵 $P$，使得 $P^{-1}AP$ 是Jordan标准型。Jordan标准型是一个块对角矩阵，其中每个块称为Jordan块。
2. Jordan块：Jordan块是一个形式为 $\lambda I + N$ 的矩阵，其中 $\lambda$ 是矩阵的特征值，$I$ 是单位矩阵，$N$ 是上三角矩阵，其对角线元素为0或者都是1。Jordan块的大小取决于特征值的代数重数和几何重数。
3. 特征值和特征向量：Jordan标准型将矩阵的特征值和特征向量整合在一起，使得特征向量更易于计算和理解。每个Jordan块对应一个特征值，而特征向量则可以从Jordan块中推导出来。

• 定理：对任意的向量 $d$， 迭代法收敛的充分必要条件是 $\rho(G)<1$

对迭代形式 $x^{(k)}=Gx^{(k-1)}+d$ 和最终解 $x^*=Gx^*+d$ 做差并递推：

$$x^{(k)}-x^*=G\left(x^{(k-1)}-x^*\right)=\cdots=G^{k}\left(x^{(0)}-x^*\right)$$

若收敛，则 $x^{(k)}-x^*\to 0$ ，即 $G^{k}\left(x^{(0)}-x^*\right)\to 0$，即 $G^{k}\to 0,k\to \infty$，由上引理得 $\rho(G)<1$

同理可证充分性

收敛性判断

用迭代矩阵的谱半径判断迭代公式是否收敛往往不容易，可以使用如下定理

如果矩阵 $G$ 的某种范数$\|G\|<1$，则

• 方程组 $x=Gx+d$ 的解 $x^*$ 存在且唯一；

• 对于迭代公式 $x^{(k+1)}=Gx^{(k)}+d$ 有

  $$\lim\limits_{k\to\infty}x^{(k)}=x^*,\quad\forall x^{(0)}\in\mathbb{R}^n$$

并且下列两式成立

$$\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|^k}{1-\|G\|}\left\|x^{(1)}-x^{(0)}\right\| \tag{2.3.6}$$

$$\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|}{1-\|G\|}\left\|x^{(k)}-x^{(k-1)}\right\|\tag{2.3.7}$$

对于第一点，因 $\|G\| < 1$ , 根据定理算子范数的估计，可知矩阵 $I-G$非奇异， 其中 $I$ 是单位矩阵

对于第二点，由 $\|G\| < 1$ ，则 $\rho(G)<1$ ，即迭代法收敛，则 $\lim\limits_{k\to\infty}x^{(k)}=x^*$

则对 $x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)$ 两边取范数，由$\|x+y\|\leq\|x\|+\|y\|$ 得：

$$\begin{aligned}\left\|x^{(k)}-x^*\right\|&\leqslant\sum_{i=k}^\infty\left\|x^{(i)}-x^{(i+1)}\right\|=\sum_{i=k}^\infty\left\|G^i\left(x^{(0)}-x^{(1)}\right)\right\|\\&\leqslant\sum_{i=k}^\infty\|G\|^i\left\|x^{(0)}-x^{(1)}\right\|\\&=\frac{\|G\|^k}{1-\|G\|}\left\|x^{(1)}-x^{(0)}\right\|,\end{aligned}$$

PS. $x^{(k)}-x^*=\sum_{i=k}^\infty\left(x^{(i)}-x^{(i+1)}\right)$ 这个等式上面证过

对于2.3.7

$$x^{(k)}-x^*=G\left(x^{(k-1)}-x^*\right)=G\left(x^{(k-1)}-x^{(k)}\right)+G\left(x^{(k)}-x^*\right).$$

则：

$$\begin{aligned}\left\|x^{(k)}-x^*\right\|&\leqslant\|G\|\left\|x^{(k-1)}-x^{(k)}\right\|+\|G\|\left\|x^{(k)}-x^*\right\|,\\&\left\|x^{(k)}-x^*\right\|\leqslant\frac{\|G\|}{1-\|G\|}\left\|x^{(k)}-x^{(k-1)}\right\|.\end{aligned}$$

Jacobi 迭代法

若系数矩阵 $A=(a_{ij})\in\mathbb{R}^{n\times n}$ 满足条件 $a_{ii}\neq0(i=1,2,\cdots,n)$，则将 $A$ 分解为 $A=D+L+U$

其中 $D$ 为对角矩阵，$L$ 为严格下三角矩阵，$U$ 为严格上三角矩阵

相当于 $D=diag(a_{ii})$ 且显然 $D^{-1}$ 存在

对于一般形式，取 $N=D,P=-(L+U)$ ，得迭代式：

$$x^{(k+1)}=-D^{-1}(L+U)x^{(k)}+D^{-1}b\quad(k=0,1,\cdots)$$

其中 $x^{(0)}\in\mathbb{R}^n$ 任取，迭代矩阵：

$$G_J=-D^{-1}(L+U)$$

此外，$d=D^{-1}b$

由 $D^{-1}=diag(a_{ii}^{-1})$，分量形式为：

$$x_i^{(k+1)}=(-\sum_{j=1}^na_{ij}x_j^{(k)}+b_i)/a_{ii}\quad(i=1,2,\ldots,n;k=0,1,\ldots)\\ x_i^{(k+1)}=-\sum_{j=1}^n\frac{a_{ij}}{a_{ii}}x_j^{(k)}+\frac{b_i}{a_{ii}}\quad(i=1,2,\ldots,n;k=0,1,\ldots)$$

• 由定理：对任意的向量 $d$， 迭代法收敛的充分必要条件是 $\rho(G)<1$，得： Jacobi 迭代法收敛的充分必要条件是 $\rho(G_J)<1$

• 由定理：矩阵 $G$ 的某种范数$\|G\|<1$，则方程组 $x=Gx+d$ 的解 $x^*$ 存在且唯一，得 ：如果 $\|G_{J}\|<1$, 则 Jacobi 迭代法收敛.

引理：

若矩阵 $A\in\mathbb{R}^{n\times n}$ 是主对角线按行 (或按列) 严格占优阵，则 $A$ 是非奇异矩阵.

主对角线按行 (或按列) 严格占优阵，即满足：

$$|a_{ii}|>\sum_{j=1\atop j\neq i}^n|a_{ij}|\:,\quad i=1,2,\cdots,n.$$

即主对角线元素大于该行非对角线元素绝对值之和

如果系数矩阵$A$是主对角线按行(或按列) 严格占优阵，则用Jacobi迭代法求解必收敛.

若 $A$ 是对称正定矩阵，则雅可比迭代收敛的充要条件是 $2D - A$ 也是对称正定矩阵.

Gauss-Seidel迭代法

在迭代法一般形式中，取$N=D+L,P=-U$，形成以下的迭代公式

$$x^{(k+1)}=-(D+L)^{-1}Ux^{(k)}+(D+L)^{-1}b\quad(k=0,1,\cdotp\cdotp\cdotp)$$

其中：

$$G_G=-(D+L)^{-1}U\\ d=(D+L)^{-1}b$$

使用时为了避免求 $D+L$ 的逆矩阵，可以通过这种方式改写：

$$\begin{aligned}&(D+L)x^{(k+1)}=-Ux^{(k)}+b\\&Dx^{(k+1)}=-Lx^{(k+1)}-Ux^{(k)}+b\\&x^{(k+1)}=-D^{-1}Lx^{(k+1)}-D^{-1}Ux^{(k)}+D^{-1}b\end{aligned}$$

分量形式：

$$x_i^{(k+1)}=-\sum_{j=1}^{i-1}\frac{a_{ij}}{a_{ii}}x_j^{(k+1)}-\sum_{j=i+1}^n\frac{a_{ij}}{a_{ii}}x_j^{(k)}+\frac{b_i}{a_{ii}}\quad(i=1,2,\ldots,n;k=0,1,\ldots)$$

• 收敛的充分必要条件是 $\rho(G_G)<1$

• $\|G_G\|<1$ 则GS法收敛

• 系数矩阵 $A$ 主对角线严格占优，则GS法必收敛

• 系数矩阵 $A$ 正定，则GS法必收敛（正定则所有特征值非负）

Gauss-Seidel迭代法相对于Jacobi迭代法的优点

在Jacobi迭代过程中，在计算解向量 $x^{(k+1)}$ 的第 $i$ 分量 $x_i^{(k+1)}$ 时，前面的 $i-1$ 个分量的新值 $x_1^{(k+1)},\quad x_2^{(k+1)},\cdots,x_{i-1}^{(k+1)}$ 已经计算出来,从而在计算 $x_i^{(k+1)}$ 时，前面的分量用新值 $x_1^{(k+1)},\quad x_2^{(k+1)},\cdots,x_{i-1}^{(k+1)}$ 后面的分量还没有计算出来，利用旧值 $x_{i+1}^{(k)},\quad x_{i+2}^{(k)},\cdots,x_{n}^{(k)}$ 来计算可以节约存储，加快速度

注意观察两种方法迭代式的特点：

• Jacobi： $D\,x^{(k+1)}=\left(b-(L+U)x^{(k)}\right)$

• Gauss-Seidel：$(D+L)\,x^{(k+1)}=b-Ux^{(k)}$

我们可以看到 Jacobi 只利用了对角信息( $A\approx D$ )，表现为每轮迭代只利用上一轮的结果；

Gauss-Seidel 利用了对角+下三角的信息( $A\approx D+L$ )，表现为及时利用本轮已经更新的结果。

矩阵特征值与特征向量计算

幂法

幂法主要用于计算矩阵的按模为最大的特征值和相应的特征向量。即：

设 $n\times n$ 实矩阵 $A$ 具有 $n$ 个线性无关的特征向量$x_1,x_2,\cdotp\cdotp\cdotp,x_n$，其相应的特征值$\lambda_1,\lambda_2,\cdotp\cdotp\cdotp,\lambda_n$满足不等式

$$\mid\lambda_1\mid>\mid\lambda_2\mid\geqslant\mid\lambda_3\mid\geqslant\cdotp\cdotp\cdotp\geqslant\mid\lambda_n\mid \tag{3.1}$$

其中$A\boldsymbol{x}_i=\lambda_i\boldsymbol{x}_i(i=1,2,\cdotp\cdotp\cdotp,n)$。现在要求出 $\lambda_1$ 和相应的特征向量。

任取一 $n$ 维非零向量$u_0$，从 $u_0$ 出发，按照如下的递推公式

$$u_k=Au_{k-1}\quad(k=1,2,\cdotp\cdotp\cdotp)\tag{3.2}$$

可产生一个向量序列 $\{u_k\}$，分析这一序列的收敛情况，可从中找出计算 $\lambda_1$ 和相应的特征向量的方法。

因 $n$ 维向量组 $x_1,x_2,\cdotp\cdotp\cdotp,x_n$ 线性无关，故对于向量 $u_0$ ，必存在唯一的不全为零的数组$\alpha_1,\alpha_2,\cdotp\cdotp\cdotp,\alpha_n$，使得

$$u_0=\alpha_1x_1+\alpha_2x_2+\cdotp\cdotp\cdotp+\alpha_nx_n=\sum_{i=1}^n\alpha_ix_i$$

所以由$A\boldsymbol{x}_i=\lambda_i\boldsymbol{x}_i(i=1,2,\cdotp\cdotp\cdotp,n)$ ， $A^ku_0=A^k\sum_1^na_ix_i=\sum_1^na_i\lambda_i^kx_i$

下面会用到

由式(3.2)可得

$$u_{k}=Au_{k-1}=A^2u_{k-2}=\cdots=A^ku_0=\\\alpha_1A^kx_1+\alpha_2A^kx_2+\cdots+\alpha_nA^kx_n=\\\alpha_1\lambda_1^kx_1+\alpha_2\lambda_2^kx_2+\cdots+\alpha_n\lambda_n^kx_n=\\\lambda_1^k\left[\alpha_1x_1+\alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^kx_2+\cdots+\alpha_n\left(\frac{\lambda_n}{\lambda_1}\right)^kx_n\right]\tag{3.3}$$

设 $\alpha_1\neq0$。由于有式(3.1)成立，故从式(3.3)可看出，当 $k$ 充分大时，有

$$u_k\approx\lambda_1^k\alpha_1x_1$$

$$u_k\approx\lambda_1^k\alpha_1x_1=\lambda_1(\lambda_1^{k-1}\alpha_1x_1)$$

所以实际上 $\lambda_1^{k-1}\alpha_1x_1$ 也可以看作一个特征向量

因为矩阵的特征向量与任何一个非零常数相乘之后仍是该矩阵的属于同一个特征值的特征向量，所以，当 $k$ 充分大时，由迭代公式(3.2)产生的 $u_k$ 可近似地作为矩阵 $A$ 的属于 $\lambda_1$ 的特征向量。迭代公式(3.2)实质上是

$$u_k=A^ku_0$$

归一化

实际计算时，为了避免迭代向量 $u_k$ 的模过大(当|$\lambda_1|>1$)或过小(当$|\lambda_1|<1$),通常每迭代一次都对 $u_k$ 进行归一化，使其范数 ($\|\cdot\|_2$或$\|\cdot\|_\infty$) 等于1。因此，实际使用的迭代公式是

$$\begin{cases}y_{k-1}=\frac{u_{k-1}}{\parallel u_{k-1}\parallel}\\u_k=Ay_{k-1}\end{cases}(k=1,2,\cdotp\cdotp\cdotp)\tag{3.4}$$

范数为1，即：

$$\|y_k\|=\|\frac{u_k}{\|u_k\|}\|=\frac{\|u_k\|}{\|u_k\|}=1$$

由 （3.4）:

$$u_{k}=\frac{Au_{k-1}}{\parallel u_{k-1}\parallel}=\frac{A^{2}u_{k-2}}{\parallel Au_{k-2}\parallel}=\cdotp\cdotp\cdotp=\frac{A^{k}u_{0}}{\parallel A^{k-1}u_{0}\parallel}$$

又：$A^ku_0=\sum_1^na_i\lambda_i^kx_i$：

$$\begin{aligned}\mathbf{y}_{k}&=\frac{\boldsymbol{u}_k}{\parallel\boldsymbol{u}_k\parallel}=\frac{\boldsymbol{A}^k\boldsymbol{u}_0}{\parallel\boldsymbol{A}^k\boldsymbol{u}_0\parallel}=\\&\left(\frac{\lambda_1}{\mid\lambda_1\mid}\right)^k\frac{\alpha_1\boldsymbol{x}_1+\alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^k\boldsymbol{x}_2+\cdots+\alpha_n\left(\frac{\lambda_n}{\lambda_1}\right)^k\boldsymbol{x}_n}{\left\|\alpha_1\boldsymbol{x}_1+\alpha_2\left(\frac{\lambda_2}{\lambda_1}\right)^k\boldsymbol{x}_2+\cdots+\alpha_n\left(\frac{\lambda_n}{\lambda_1}\right)^k\boldsymbol{x}_n\right\|}\end{aligned}\tag{3.5}$$

注意这个上下提出 $\lambda_1^k$，范数中提出后要加绝对值

由于 $\mid\lambda_1\mid>\mid\lambda_2\mid\geqslant\mid\lambda_3\mid\geqslant\cdotp\cdotp\cdotp\geqslant\mid\lambda_n\mid$，故当 $k\to \infty$ :

• $\lambda_1>0$，则 $y_k\to\frac{\alpha_1x_1}{\parallel\alpha_1x_1\parallel}$ ，即 $y_k$ 有确定的极限
• $\lambda_1<0$，则 $y_k$ 的各个分量之模有确定的极限，其各个分量的正负号每迭代一次就改变一次

$\lambda_1$ 计算方法1：2范数

关于 $\lambda _1$ 的 计 算，有两种方法可供选择。第 一 种方法，在式(3.4) 中使用范数 $\| \cdot \| _2$ ，并且令

$$\beta_k=y_{k-1}^\mathrm{T}u_k\tag{3.6}$$

那么，由于 $u_k=Ay_{k-1}$ ，并根据式(3.5)可得

$$\lim_{k\to\infty}\beta_k=\frac{\alpha_1x_1^\mathrm T}{\parallel\alpha_1x_1\parallel_2}A\frac{\alpha_1x_1}{\parallel\alpha_1x_1\parallel_2}=\frac{\alpha_1^2x_1^\mathrm Tx_1}{\parallel\alpha_1x_1\parallel_2^2}\lambda_1=\lambda_1$$

注意这里使用的二范数的定义$\parallel x_1\parallel_2^2=x_1^\mathrm{T}x_1$, $\parallel\alpha_1x_1\parallel_2^2=\alpha_1^2x_1^\mathrm{T}x_1$

把式(3.4)和式(3.6)结合在一起，得到第一种幂法迭代格式：

$$\begin{cases}\text{任取非零向量 }u_{0}\in\mathbf{R}^{n}\\\eta_{k-1}=\sqrt{\boldsymbol{u}_{(k-1)}^{\mathrm T}\boldsymbol{u}_{k-1}}\\\boldsymbol{y}_{k-1}=\boldsymbol{u}_{k-1}/\boldsymbol{\eta}_{k-1}\\\boldsymbol{u}_{k}=\boldsymbol{A}\boldsymbol{y}_{k-1}\\\boldsymbol{\beta}_{k}=\boldsymbol{y}_{k-1}^{\mathrm{T}}\boldsymbol{u}_{*}\\(k=1,2,\cdotp\cdotp\cdotp)&&\end{cases}$$

$\eta_{k-1}=\sqrt{\boldsymbol{u}_{(k-1)}^{\mathrm T}\boldsymbol{u}_{k-1}}$ 就是2范数

$\boldsymbol{y}_{k-1}=\boldsymbol{u}_{k-1}/\boldsymbol{\eta}_{k-1}$ 对应归一化的过程

当$|\beta_k-\beta_{k-1}|/|\beta_k|\leqslant\varepsilon$ (允许误差)时，选代终止，以当前的 $β_k$ 作为 $\lambda_1$ 的近似值，以 $y_{k-1}$作为 $A$ 的属于 $\lambda_1$ 的特征向量。

不用 $|\beta_k-\beta_{k-1}|$ 作为判断依据，避免 $\beta_k$ 本身就很小的情况

$\lambda_1$ 计算方法2：无穷范数

令：

$$\beta_k=\frac{e_r^\mathrm{T}u_k}{e_r^\mathrm{T}y_{k-1}}$$

这里假定 $u_{k-1}$的第 $r$ 个分量为模最大的分量，当 $k$ 足够大之后，$r$ 保持定值；$e_r$ 是 $n$ 维基本单位向量，它的第 $r$ 个分量为1，其余分量为零。由于 $u_k=Ay_{k-1}$，并根据式(3.5)可得

$$\lim_{k\to\infty}\beta_k=\frac{e_r^\mathrm{T}Ax_1}{e_r^\mathrm{T}x_1}=\lambda_1$$

把式(3.4)和(3.8)结合在一起，得到第二种幂法迭代格式：

$$\begin{cases}&\text{任取非零向量 }u_0=(h_1^{(0)},\cdotp\cdotp\cdotp,h_n^{(0)})^\mathrm{T}\\&|h_r^{(k-1)}|=\max_{1\leqslant j\leqslant n}|h_j^{(k-1)}|\\&y_{k-1}=u_{k-1}/|h_r^{(k-1)}|\\&u_k=Ay_{k-1}=(h_1^{(k)},\cdotp\cdotp\cdotp,h_n^{(k)})^\mathrm{T}\\&\beta_k=\operatorname{sgn}(h_r^{k-1})h_r^{(k)}\\&(k=1,2,\cdots)\end{cases}\tag{3.9}$$

终止迭代的控制也用$|\beta_k-\beta_k|/|\beta_k|\leqslant\varepsilon$,当前的 $\beta_k$ 和 $y_k-1$即分别作为 $\lambda_1$ 和与其相应的特征向量。在迭代格式(3.9)中，$|h_r^{k-1}|=\|u_{k-1}\|_{\infty}$ ，$sgn(h_r^{(k-1)})=e_r^{\mathrm{T}}y_{k-1},h_r^{(k)}=e_r^{\mathrm{T}}u_k$。

收敛速度

幂法的收敛速度与比值$\left|\frac{\lambda_2}{\lambda_1}\right|$ [式(3.1)的情况] 或$\left|\frac{\lambda_{m+1}}{\lambda_1}\right|$ [式(3.10)的情况] 有关，比值越小，收敛速度越快。

此外，当矩阵$\boldsymbol{A}$没有$n$个线性无关的特征向量，幂法仍然可以使用，但收敛速度特别慢，应改用其他方法。

反幂法

幂法主要用于计算矩阵的按模为最小的特征值 $\lambda_n$ 和相应的特征向量。

因 $A$ 非奇异，故$\lambda_i\neq0(\boldsymbol{i}=1,2,\cdotp\cdotp\cdotp,n)$。由

$$Ax_i=\lambda_ix_i$$

得

$$A^{-1}x_i=\frac1{\lambda_i}x_i\quad(i=1,2,\cdotp\cdotp\cdotp,n)$$

所以，$\frac1{\lambda_n}$ 是矩阵 $A^{-1}$ 的按模最大的特征值，$x_n$ 是 $A^{-1}$ 的属于$\frac1{\lambda_n}$的特征向量。于是，对矩阵 $A^{-1}$ 使用幂法迭代格式(3.7)或式(3.9)就可求出 $\frac1{\lambda_n}$ (从而求出 $\lambda_n$) 以及相应的特征向量。

平移法求任意特征值及其加速

设已知数 $\mu$ 是 $n\times n$ 矩阵 $A$ 的某个特征值 $\lambda_i$的近似值，并满足

$$0< \mid \lambda , - \mu \mid < \mid \lambda _i- \mu \mid, 1\leqslant i\leqslant n, i\neq s$$

于是，对矩阵$A-\mu I$ 实行反幂法迭代，就可求出 $A$ 的特征值 $\lambda_i$ 。及其相应的特征向量。其中 $\mu$ 称为平移量。具体说，就是把反幂法迭代格式中的 $Au_k=y_{k-1}$改为

$$(A-\mu I)u_k=y_{k-1}$$

当 $k$ 足够大时，取 $\lambda_s\approx\beta_k^{-1}+\mu$ ，其相应的特征向量近似地等于$y_{k-1}$。

由于幂法和反幂法的迭代是否收敛依赖于特征值的分布情况，因此实际使用时很不方便， 特别是不适合于自动计算。

这里我们使用反幂法迭代而非幂法，因为我们不容易通过平移来确保平移后特征值的绝对值最大，但容易使其尽量小

平移加速收敛

为方便说明，假设 $A$ 的特征值均为实数且大于零，则$B=A-\mu I$，平移前的收敛速度 $\left | \frac {\lambda _2}{\lambda _1}\right |$，平 移后的收敛速度 $\left | \frac {\lambda _2( B) }{\lambda _1( B) }\right |$。 对任意的$0<\mu<\lambda_n$，我们有：

$$|\frac{\lambda_2(B)}{\lambda_1(B)}|<\left|\frac{\lambda_2}{\lambda_1}\right|$$

也就意味着平移可以加速收敛

Rayleigh( 瑞利 ）商加速

Reference

https://blog.csdn.net/HG0724/article/details/122048108

https://zhuanlan.zhihu.com/p/389389672