2022-2023 春学期矩阵与数值分析 C4 逐次逼近法

原文

C4 逐次逼近法

4.1 解线性方程组的迭代法

简单迭代法迭代格式

可将线性方程组变形

\[Ax=b \Leftrightarrow x=Bx+f \]

其中 B 为迭代矩阵，且

\[B\in R^{n\times n},f\in R^n,x\in R^n \]

迭代法：称使用

\[x^{(k+1)}=Bx^{k}+f\;(k=0,1,2,\cdots) \]

求解的方法为迭代法，也称迭代过程或迭代格式.

迭代法收敛、发散：如果对任意 \(x^{(0)}\) ，都有 \(\lim\limits_{k\to\infty}x^{(k)}=x^*\)。即

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

称该迭代法收敛，否则称迭代法发散。

Jacobi，Gauss-seidel 迭代法等

简单迭代法-方程组变形

线性方程组：

\[\left\{ \begin{array}{c} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1\\ \vdots\\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n \end{array} \right. \]

等价方程组：\(a_{ii}\neq0\)

\[\left\{ \begin{array}{l} x_1=\frac{1}{a_{11}}(b_1-a_{12}x_2-\cdots-a_{1n}x_n)\\ \qquad\qquad\qquad\qquad\qquad\vdots\\ x_i=\frac{1}{a_{ii}}(b_i-a_{i1}x_1-\cdots-a_{ii-1}x_{i-1}-a_{ii+1}x_{i+1}-\cdots-a_{in}x_n)\\ \qquad\qquad\qquad\qquad\qquad\vdots\\ x_n=\frac{1}{a_{nn}}(b_a-a_{n1}x_1-\cdots-a_{nn-1}x_{n-1}) \end{array} \right. \]

Jacobi 迭代法：

\[x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum\limits^{i-1}_{j=1}a_{ij}x_j^{(k)}-\sum\limits^n_{j=i+1}a_{ij}x_j^{(k)})\;(i=1,2,\cdots,n) \]

在 Jacobi 迭代过程中，对已经算出来的信息未加以充分利用，比如在计算 \(x^{k+1}_2\) 时 \(x^{k+1}_1\) 已经算出，可是使用 Seidel 技巧加以改进，得到 Gauss-Seidel 迭代法：

\[x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum\limits^{i-1}_{j=1}a_{ij}x_j^{(k+1)}-\sum\limits^n_{j=i+1}a_{ij}x_j^{(k)})\;(i=1,2,\cdots,n) \]

矩阵分裂：可以将系数矩阵 A 分裂为 A=D-L-U 的形式，具体如下所示：

\[A=D-L-U\\ A= \begin{pmatrix} a_{1,1}&a_{1,2}&\cdots&a_{i,j}&\cdots&a_{1,n}\\ a_{2,1}&a_{2,2}&\cdots&\cdots&\cdots&a_{2,n}\\ \vdots&\ddots&\ddots&\vdots&\vdots&\vdots \\ a_{i,1}&\cdots&a_{i,j-1}&a_{i,j}&\cdots&a_{i,n} \\ \vdots&&\vdots&\ddots&\ddots&\vdots \\ a_{n,1}&\cdots&a_{n,j-1}&\cdots&a_{n,n-1}&a_{n,n} \end{pmatrix}\\ D=\mathrm{diag(a_{1,1},a_{2,2},\cdots,a_{n,n})}\\ L= \begin{pmatrix} 0&&&&&0\\ -a_{2,1}&0\\ \vdots&\ddots&\ddots\\ -a_{j,1}&\cdots&-a_{j,j-1}&0\\ \vdots&&\vdots&\ddots&\ddots\\ -a_{n,1}&\cdots&-a_{n,j-1}&\cdots&-a_{n,n-1}&0 \end{pmatrix}\\ U= \begin{pmatrix} 0&-a_{1,2}&\cdots&-a_{1,j}&\cdots&-a_{1,n}\\ &\ddots&\ddots&\vdots&&\vdots\\ &&0&-a_{j-1,j}&\cdots&-a_{j-1,n}\\ &&&\ddots&\ddots&\vdots\\ &&&&0&-a_{n-1,n}\\ 0&&&&&0\\ \end{pmatrix}\\ \]

Jacobi、Gauss-Seidel 迭代法矩阵格式

由 A=D-L-U，得 Dx=(L+U)x+b 从而有：

\[x=D^{-1}(L+U)x+D^{-1}b \]

则 Jacobi 迭代法可写成：

\[x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b;\;k=0,1,2,\cdots \]

由 A=D-L-U，得 (D-L)x=Ux+b 从而

\[x=(D-L)^{-1}Ux+(D-L)^{-1}b \]

则 Gauss-Seidel 迭代法可以写成：

\[x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b;\;k=0,1,2,\cdots \]

并不是对任何情况，G-S 迭代都比 Jacobi 迭代收敛速度快

例题

*迭代改善法（应该不会考）

迭代法的收敛性

迭代法收敛性定理（谱半径）：迭代法 \(x^{(k+1)}=Bx^{(k)}+f\) 对任意 \(x^{(0)}\) 和 \(f\) 均收敛的充要条件为：

\[\rho(B)<1 \]

迭代法收敛性定理（范数）：若迭代矩阵任意范数 \(||B||<1\)，则迭代法收敛，且有

\[||x^{(k)}-x^*||\leq\frac{||B||^k}{1-||B||}||x^{(1)}-x^{(0)}|| \]

定理： \(\lim\limits_{k\to\infty}\epsilon^{(k)}=0\)（即迭代收敛 \(x_i^{(k)}\to x_i^*,i=1,2,\cdots,n\)）的充分必要条件是 \(\lim\limits_{k\to\infty}B^k=0_{n\times n}\Leftrightarrow \rho(B)<1\)

迭代法收敛性-特殊线性方程组

针对某些特殊方程组，从方程组本身就可以判定其收敛性，不必求迭代矩阵的特征值或范数

严格（不可约）对角占优矩阵
- 定义：如果矩阵 \(A=(a_{ij})_{n\times n}\) 的元素满足不等式 \(|a_{ii}|\geq\sum\limits_{j=1,j!=i}^n|a_{ij}|,\;(i=1,2,\cdots,n)\)，则称 A 为对角占优阵；如果所有严格不等式（即不等式不取等号）均成立，称 A 为严格对角占优矩阵。
- 定理：严格对角占优阵 A 为非奇异矩阵，即 \(\det(A)\neq0\)
- 定理：若线性方程组 Ax=b 中的 A 为严格对角占优矩阵，则 Jacobi 法和 Gauss-Seidel 法均收敛。
对称正定矩阵

若 A 为对称正定矩阵：
- 若 2D-A 为正定矩阵，则 Jacobi 迭代法收敛
- Gauss-Seidel 迭代法收敛

例题

4.2 非线性方程迭代解法

迭代法收敛

几何表示

收敛性定理：设迭代函数 \(\varphi(x)\) 满足

当 \(x\in[a,b]\) 时，\(a\leq\varphi(x)\leq b\)
存在正数 \(0<L<1\)，对任意 \(x\in[a,b]\) 均有 \(|\varphi'(x)|\leq L\)（难以验证）

则 \(x=\varphi(x)\) 在 \([a,b]\) 内存在唯一根 \(a\)，且对任意初始值 \(x_0\in[a,b]\) ，迭代法

\[x_{k+1}=\varphi(x_k)\;(k=0,1,2,\cdots) \]

收敛于 \(a\)，且

\[|x_k-a|\leq\frac{L}{1-L}|x_k-x_{k-1}|\leq\frac{L^k}{1-L}|x_1-x_0| \]

L 较小，相邻两次计算值的偏差 \(|x_k-x_{k-1}|\leq\delta\) （事先给定的精度），迭代过程就可以终止，\(x_k\) 就可以作为 a 的近似值

L 的大小可作出估计时，估计所需要的迭代次数

\[|x_k-a|\leq\frac{L^k}{1-L}|x_1-x_0|<\delta\rightarrow k\geq \log_L\frac{\delta(1-L)}{|x_1-x_0|} \]

使用迭代法时往往在根 a 的附近进行，假定 \(\varphi'(x)\) 在 a 的附近连续且满足 \(|\varphi'(a)|<1\)，则一定存在 a 的某个领域 \(S:|x-a|\leq\delta\) ，\(\varphi(x)\) 在 S 上满足收敛性定理中的条件，所以在 S 中任取初始值 \(x_0\)，迭代格式

\[x_k=\varphi(x_{k-1}) \]

收敛于方程的根 a，即 \(f(a)\equiv0\)，称此收敛为局部收敛。

误差估计式：

\[|x_k-a|\leq\frac{L^k}{1-L}|x_1-x_0| \]
中 L 或 \(|\varphi'(x)|\) 在 \([a,b]\) 上的值越小，迭代的收敛速度就越快。L < 1 且接近于 1 时，迭代法虽然收敛，但收敛速度很慢。

为了使收敛速度有定量的判断，特介绍收敛速度的阶的概念，作为判断迭代法收敛速度的重要标准。

设迭代格式 \(x_{k+1}=\varphi(x_k)\)，当 \(k\to\infty\) 时，\(x_{k+1}\to a\) 记 \(e_k=x_k-a\)

迭代法收敛阶：若存在实数 \(p\geq1\) 和 \(c>1\) 满足

\[\lim\limits_{k\to\infty}\frac{|e_{k+1}|}{|e_k|^p}=c \]

则称迭代法是 p 阶收敛，当 p=1 时，称线性收敛，当 p>1 时称超线性收敛，当 p=2 时称平方收敛。

显然：

收敛阶越大迭代法的收敛速度也越快
收敛阶难以直接确定，课采用一些其他方法来确定收敛阶（Taylor 展开式）

收敛阶判定定理

p 阶收敛：如果 \(x=\varphi(x)\) 中的迭代函数 \(\varphi(x)\) 在根 a 附近满足：

\(\varphi(x)\) 存在 p 阶导数且连续
\[ \varphi'(a)=\varphi''(a)=\cdots=\varphi^{(p-1)}(a)=0,\;\varphi^{(p)}(a)\neq0 \]

则迭代法 \(x_{k+1}=\varphi(x_k)\) 是 p 阶收敛的

例题

牛顿（Newton）法：

Newton 迭代法

Newton 迭代格式：设方程 \(f(x)=0\) 的根为 \(a\)，且 \(f'(a)\neq0\)，则迭代法

\[x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}\;(k=0,1,2,\cdots) \]

至少是平方收敛，并称该迭代格式为 Newton 迭代法。

类似的，可以使用导数的近似式 \(f'(x)\approx\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}\) 替代 \(f'(x)\)，可将其代入 Newton 迭代格式，得到割线法（弦截法）。

\[x_{k+1}=x_k-\frac{(x_k-x_{k-1})f(x_k)}{f(x_k)-f(x_{k-1})} \]

割线法是超线性收敛的（\(p=\frac{1+\sqrt{5}}{2}\approx1.618\)），低于 Newton 法

此外，还有单步弦截法

\[x_{k+1}=x_k-\frac{(x_k-x_{0})f(x_k)}{f(x_k)-f(x_{0})} \]

平行线法

\[x_{k+1}=x_k-\frac{f(x_k)}{f'(x_0)} \]

例题

针对重根改进的 Newton 法：

设 \(a\) 是 \(f(x)=0\) 的 \(m\) 重根（满足 \(f(a)=f'(a)=f''(a)=\cdots=f^{(m-1)}(a)=0,\;f^{(m)}\neq0\) ），其中 \(m\geq2\) 是整数，则有 \(f(x)=(x-a)^mg(x)\)，且 \(g(a)\neq0\)

此时，直接使用 Newton 法是线性收敛而不是平方收敛的

\[\varphi(x)=x-\frac{f(x)}{f'(x)}=x-\frac{(x-a)^mg(x)}{m(x-a)^{(m-1)}g(x)+(x-a)^mg'(x)} \\ \varphi'(x)=1-\frac{g(x)}{mg(x)+(x-a)g'(x)}-(x-a) [\frac{g(x)}{m\cdot g(x)+(x-a)g'(x)}]' \\ \varphi'(a)=1-\frac{g(a)}{m\cdot g(a)+0}-0=1-\frac{1}{m}\neq0 \]

为了提高收敛的阶，可取 \(\varphi(x)=x-m\frac{f(x)}{f'(x)}\)，从而 \(\varphi(a)=1-\frac{m}{m}=0\)，则新迭代法至少是平方收敛的

例题

*多根区间上的逐次逼近法(应该不会考)

这里主要使用二分法

4.3 幂法

参考幂法与反幂法 - 知乎 (zhihu.com)

幂法（用于求绝对值最大的特征值）

幂法用于计算绝对值最大的特征值和对应的特征向量，对于单根和重根情形均有效。

\[v_k=A^kv_0 \]

计算可得

\[v_k=\alpha_1\lambda_1^kx_1,\quad \alpha_1\neq0,\quad \frac{v_k}{v_{k-1}}=\lambda_1 \]

实际使用上，对幂法进行了改进，添加了规范化操作

算法过程

输入矩阵 \(A=(a_{ij})\)，初始向量 \(v_0\)，允许的误差 \(\epsilon\)，最大迭代次数 \(N\)
置 \(k=1,\quad \mu=0\)
求 \(\max\{v_k\},\quad\max\{v_k\}\Rightarrow a\)
计算 \(u=v/a,\quad v=Au\)，置 \(\max\{v_k\}=\lambda\)
若 \(|\lambda-\mu|<\epsilon\)，输出 \(\lambda,u\)，停机，否则跳转至 6
若 \(k<N\)，置 \(k+1\Rightarrow k,\lambda\Rightarrow\mu\)，转到 3，否则输出失败信息，停机

例题

利用幂法计算矩阵 A 的主（最大）特征值和相应的特征向量，其中 \(A=\begin{pmatrix}2&3&2\\10&3&4\\3&6&1\end{pmatrix}\)，取 \(v_0=(0,0,1)^T\)，迭代2次

解：

取 \(u_0=v_0=(0,0,1)^T\)，则 \(u_1=Av_0=(2,4,1)^T\)，且 \(\max\{u_1\}=-4\)

即 \(a=4\)，于是 \(v_1=u_1/a=(0.5,1,0.25)^T\)

\(u_2=Av_1=(4.5,9,7.75)^T\) 且 \(\max\{u_2\}=9\)

即 \(a=9\)，于是 \(v_2=u_2/a\approx(0.5,1,0.8611)\)，可得 \(\lambda_1\approx9,x_1\approx(0.5,1.0,0.8611)^T\)

注：经过 8 次迭代可得 \(\lambda_1=11,x_1=(0.5,1.0,0.75)^T\)

反幂法

反幂法用于计算矩阵 A 按模最小的特征值及对应的特征向量

算法流程：

任取初始向量 \(v_0=u_0\neq0\)
计算 \(v_k=A^{-1}u_{k-1},\;(k=1,2,\cdots)\)
\(\max\{v_k\}\Rightarrow m_k;\quad u_k=v_k/m_k\)
如果 \(k\) 从某时以后有 \(\frac{(v_k)_j}{(v_{k-1})_j}\approx c \;\;(常数)(j=1,2,\cdots,n)\)，则取 \(\lambda_n\approx\frac{1}{c}\)；而 \(u_k\) 就是与 \(\lambda_k\) 对应的特征向量

例题

用反幂法求矩阵 \(A=\begin{pmatrix}3&2\\4&5\end{pmatrix}\) 的按模最小的特征值和特征向量，精确至 6 位小数，迭代 1 次

解：

取 \(v_0=u_0=(1,1)^T\)，由 \(v_1=A^{-1}u_0\)，即 \(Av_1=u_0\)

得 \(v_1=(0.428571,-0.142857)^T\\m_1=\max\{v_1\}=0.428571\)

归一化得 \(u_1=(1.000000,-0.333333)^T\)

所以迭代 1 次得到的按模最小的特征值为 \(1/0.428571\)，对应的特征向量为 \(u_1=(1.000000,-0.333333)^T\)

用反幂法求矩阵 \(A=\begin{pmatrix}3&2\\4&5\end{pmatrix}\) 的过程中，得到 \(u^{(8)}=\begin{pmatrix}1.000000\\-0.999999\end{pmatrix}\)，则 A 的模最小的特征值的近似值为：

解：

由于 \(c=\max(u^{(8)}=1)\)，则按模最小的特征值的近似值为 \(\lambda_1=\frac{1}{c}=1.000000\)

或者

依据特征值的概念与性质，可知有 \(Ax=\lambda x\)，则

\[Au^{(8)}=\lambda_nu^{(8)}\approx \begin{pmatrix} 1\\-1 \end{pmatrix} \]

则，可认为按模最小的特征值的近似值为 1。

4.4 共轭梯度法

共轭梯度法是以共轭方向作为搜索方向的一类算法，最初用于求解线性方程组，后来用于求解无约束最优化问题，是一种重要的数学优化方法。

*最速下降法（应该不会考）

共轭梯度法

下降方向 p 为上步下降方向与当前残量的线性组合，与 \(\beta\) 有关，步长为 \(\alpha\)

共轭梯度法最多迭代 n 步即可求得问题的精确解

算法的流程：

对于所给的初始向量 \(x_0\)，\(r_0=p_0=b-Ax_0\)

\[\alpha_k=(r_k,r_k)/(Ap_k,p_k)\\ x_{k+1}=x_k+\alpha_kp_k\\ r_{k+1}=r_k-\alpha_kAp_k\\ \beta_k=(r_{k+1},r_{k+1})/(r_k,r_k)\\ p_{k+1}=r_{k+1}+\beta_kp_k\\ (k=0,1,2,\cdots) \]

\((r_k,r_k)\) 表示内积运算。