数值方法笔记3：线性和非线性方程组求解

前置知识1：矩阵范数

矩阵范数的性质：

(1) 若$A\ne 0$ , 那么$\Vert A\Vert >0$; 若 $A=0$ , 那么$\Vert A\Vert =0$ .
(2) 对于$\alpha \in \mathrm{R}$,$\Vert \alpha A\Vert =|\alpha |\Vert A\Vert$ .
(3) $\Vert A+B\Vert \le \Vert A\Vert +\Vert B\Vert$.
(4) $\Vert AB\Vert \le \Vert A\Vert \Vert B\Vert$ .

常见的矩阵范数：

$$\begin{array}{rl}& \Vert A\Vert =\underset{\Vert x{\Vert }_1}{\max }\Vert Ax\Vert \\ & \Vert A{\Vert }_1=\underset{\Vert x{\Vert }_1=1}{\max }\Vert Ax{\Vert }_1=\underset{k}{\max }\sum _{i=1}^n\left|a_{ik}\right|(\text{列和范数})\\ & \Vert A{\Vert }_2=\underset{\Vert x{\Vert }_2=1}{\max }\Vert Ax{\Vert }_2=\sqrt{{\lambda }_1},{\lambda }_1\text{ 是}{A}^{T}A\text{的最大特征值（}A\text{的最大奇异值）}.\\ & \Vert A{\Vert }_{\infty }=\underset{\mid x{\Vert }_{\infty }=1}{\max }\Vert Ax{\Vert }_{\infty }=\underset{i}{\max }\sum _{k=1}^n\left|a_{ik}\right|(\text{行和范数})\\ & \Vert A{\Vert }_F={\left(\sum _{i=1,k=1}^n{\left|a_{ik}\right|}^2\right)}^{1/2}\end{array}$$

而对于向量范数：

$$\Vert x{\Vert }_p={\left(\sum _{k=1}^n{\left|x_k\right|}^p\right)}^{1/p}$$

范数$\Vert \ast {\Vert }_a$和范数$\Vert \ast {\Vert }_b$等价：
$$\begin{array}{r}{c}_1\Vert A{\Vert }_a\le \Vert A{\Vert }_b\le c_2\Vert A{\Vert }_a\\ c_1^{\prime }\Vert A{\Vert }_b\le \Vert A{\Vert }_a\le c_2\Vert A{\Vert }_b\end{array}$$

当$A$的矩阵范数$\Vert A\Vert <1$，则$I\pm A$是非奇异可逆矩阵：
$$\Vert (I\pm A{)}^{-1}\Vert \le \frac{1}{1-\Vert A\Vert }$$

前置知识2：舒尔补

$$\begin{array}{rl}A& =\left[\begin{array}{cc}B& C\\ D& E\end{array}\right],\det (B)\ne 0\\ A& =\left[\begin{array}{cc}I& 0\\ D{B}^{-1}& I\end{array}\right]\left[\begin{array}{cc}B& 0\\ 0& E-D{B}^{-1}C\end{array}\right]\left[\begin{array}{cc}I& B^{-1}C\\ 0& I\end{array}\right]\end{array}$$

前置知识3：可约矩阵

定义：
如果通过行列变换可以变成这种形式：
$$\mathrm{PAQ}=\left[\begin{array}{cc}\mathbf{F}& 0\\ \mathbf{G}& \mathbf{H}\end{array}\right]或\left[\begin{array}{cc}\mathbf{F}& \mathbf{G}\\ 0& \mathbf{H}\end{array}\right]$$
左下角或右上角的$0$是零矩阵，则$A$是可约矩阵。

可约矩阵：
$$\left[\begin{array}{llll}2& 0& 1& 0\\ 8& 6& 7& 5\\ 4& 2& 3& 1\\ 4& 0& 3& 0\end{array}\right]\Rightarrow \left[\begin{array}{cccc}1& 2& 3& 4\\ 5& 6& 7& 8\\ 0& 0& 1& 2\\ 0& 0& 3& 4\end{array}\right]$$

不可约矩阵：
$$C_1=\left[\begin{array}{rrrrr}2& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 2\end{array}\right]$$

注意上面这个矩阵对角占优，但不是严格对角占优。

前置知识4：谱半径

$|A-\lambda I|=0$
定义：设${\lambda }_i$是$A$的特征值，
$$\rho (A)=\underset{1⩽i⩽n}{\max }\left|{\lambda }_i\right|$$
称为矩阵A的谱半径。
定理：矩阵谱半径和矩阵范数有如下关系：
(1)若A是一般方阵，则P(A)不能作为矩阵的范数；
(2)若A是一般方阵，则谱半径不超过任意一种矩阵范数，即
$$\rho (A)\le \Vert A\Vert$$
(3)若A为实对称矩阵，则谱半径可作矩阵的模，此时有$\rho (A)=\Vert A{\Vert }_2$。

1.【线性方程组】直接求解：高斯消元法($LU$分解)、$LDV$分解、$LD{L}^T$分解、$UD{U}^T$分解

1.1 高斯消元法($LU$分解)

针对线性方程组：
$$\mathrm{AX}=\mathrm{b}$$

我们可以将$\mathrm{A}$根据LU分解，分解为$\mathrm{A}=\mathrm{LU}$，其中$\mathrm{L}$和$\mathrm{U}$分别是下三角和上三角矩阵。

$$L=\left[\begin{array}{cccc}\ast & 0& 0& 0\\ \ast & \ast & 0& 0\\ \ast & \ast & \ast & 0\\ \ast & \ast & \ast & \ast \end{array}\right]U=\left[\begin{array}{cccc}\ast & \ast & \ast & \ast \\ 0& \ast & \ast & \ast \\ 0& 0& \ast & \ast \\ 0& 0& 0& \ast \end{array}\right]$$

有：

$$\begin{array}{rl}& \mathrm{AX}=\mathrm{b}\\ \Rightarrow & (\mathrm{LU})\mathrm{X}=\mathrm{b}\\ \Rightarrow & \mathrm{L}(\mathrm{UX})=\mathrm{b}\\ \Rightarrow & {\mathrm{LX}}_1=\mathrm{b}\mathrm{UX}={\mathrm{X}}_1\end{array}$$

高斯消元变换三角阵：①交换行。②行乘一个因子。③某一行加到另一行上。例子：





高斯消元法中如果碰到对角线上的元素（主元素）消元为0，需要交换行，称作pivot element。

当主元素不合适由于舍入误差可能会无法求解！！

所以要选择合适的主元素：

当对角线的元素是0，可以换主元素。

当然也可以提前换主元素。

这可以表示为：

总结为：
$$\begin{array}{rl}& Ax=bPA=LU\\ \Rightarrow & PAx=Pb=\bar{b}\\ \Rightarrow & LUx=\bar{b}\\ \Rightarrow & L\bar{x}=\bar{b}Ux=\bar{x}\end{array}$$

计算的复杂度：

乘法和除法：${\sum }_{p=1}^{N-1}(N-p)(N-p+1)=\frac{N^3-N}{3}$
减法：${\sum }_{p=1}^{N-1}(N-p)(N-p)=\frac{2N^3-3N^2+N}{6}$
这里使用了公式：$$\sum _{k=1}^{M}k=\frac{M(M+1)}{2},\sum _{k=1}^M{k}^2=\frac{M(M+1)(2M+1)}{6}$$

1.2 $LDV$分解、$LD{L}^T$分解、$UD{U}^T$分解

定理：
$$a_{kk}^{(k)}\ne 0,k=1,2,\dots ,n\iff \left|A_k\right|\ne 0,k=1,2,\dots ,n$$

其中$|A_k|$是方阵的$k$阶主子式。

进而我们可以知道$\mathrm{A}=\mathrm{LU}$,$\mathrm{L}$是下三角矩阵，$\mathrm{U}$是上三角矩阵。$\mathrm{L}$的对角线元素都是1，$\mathrm{L}$的行列式$|\mathrm{L}|$是1，所以$|\mathrm{A}|=|LU|=|\mathrm{U}|$

让$\mathrm{U}=\mathrm{DR}$，则有：

$$\mathrm{A}=\mathrm{LDR}$$

其中$\mathrm{L}$是下三角矩阵，$\mathrm{D}$是对角矩阵，$\mathrm{R}$是上三角矩阵。

当$\mathrm{A}$是对称阵，$\mathrm{A}={\mathrm{LDL}}^{\mathrm{T}}$

类似的对称阵也可以表示为$\mathrm{A}={\mathrm{UDU}}^{\mathrm{T}}$

还有就是求逆矩阵的方法：

$\left[\begin{array}{ll}A& I\end{array}\right]\stackrel{\text{ Row Transformation }}{⟶}\left[\begin{array}{ll}I& A^{-1}\end{array}\right]$

总结：直接求解线性方程组：

1.核心算法是LU分解。 $$\begin{array}{rl}PB& =LU\\ B^{-1}& =U^{-1}{L}^{-1}P\end{array}$$

2.迭代求解器可能不能收敛或计算成本较高。

1.3 误差分析（从条件数的角度）

矩阵范数回顾前置知识1.
①给$b$一个小的扰动对$x$有什么影响嘛？

$$\begin{array}{rl}& Ax=b\\ & A(x+\delta x)=b+\delta b\\ & A\delta x=\delta b\\ & \Vert \delta x\Vert =\Vert A^{-1}\delta b\Vert \le \Vert A^{-1}\Vert \delta b\Vert \\ & \Vert Ax\Vert \le \Vert A\Vert \Vert x\Vert ,\Vert x\Vert \ge \frac{\Vert Ax\Vert }{\Vert A\Vert }=\frac{\Vert b\Vert }{\Vert A\Vert }\\ & \frac{\Vert \delta x\Vert }{\Vert x\Vert }\le \Vert A\Vert \Vert A^{-1}\Vert \frac{\Vert \delta b\Vert }{\Vert b\Vert }=\mathrm{cond}(A)\frac{\Vert \delta b\Vert }{\Vert b\Vert }\end{array}$$

②给$A$一个小的扰动对$x$有什么影响嘛？

$$\begin{array}{rl}& Ax=b\\ & (A+\delta A)(x+\delta x)=b\\ & A\delta x+\delta A(x+\delta x)=0\\ & \Vert \delta x\Vert =\Vert A^{-1}\delta A(x+\delta x)\Vert \le \Vert A^{-1}\Vert \delta A\Vert \Vert (\Vert x\Vert +\Vert \delta x\Vert )\\ & \frac{\Vert \delta x\Vert }{\Vert x\Vert }\le \Vert A^{-1}\Vert \Vert \delta A\Vert \left(1+\frac{\Vert \delta x\Vert }{\Vert x\Vert }\right)\\ & \frac{\Vert \delta x\Vert }{\Vert x\Vert }\le \frac{\Vert A^{-1}\Vert \Vert \delta A\Vert }{1-\Vert A^{-1}\Vert \Vert \delta A\Vert }=\frac{\Vert A^{-1}\Vert \Vert A\Vert \frac{\Vert \delta A\Vert }{\Vert A\Vert }}{1-\Vert A^{-1}\Vert \Vert A\Vert \frac{\Vert \delta A\Vert }{\Vert A\Vert }}=\frac{cond(A)|\frac{\Vert \delta A\Vert }{\Vert A\Vert }}{1-cond(A)|\frac{\Vert \delta A\Vert }{\Vert A\Vert }}\end{array}$$

这里需要假设$1-cond(A)\frac{\Vert \delta A\Vert }{\Vert A\Vert }\ge 0$或$\delta A$足够小。其中$cond(A)=\Vert A^{-1}\Vert \Vert A\Vert$称为条件数。

当条件数很大时矩阵是病态的。例如:


其他判断条件数的方法

1. 当两行中的对应元素的比率非常接近时，cond(A)可能很大。
2. 元素之间的差异很大，cond(A)也可能很大。
3. 对A或b做一个小的扰动，然后解方程。如果解差很大，矩阵就没有条件。

于是我们就想把病态的矩阵转为非病态的矩阵：左乘一个矩阵，改变稳定性。

找到合适的${\tilde{A}}^{-1}$是关键

2. 【线性方程组】间接迭代求解：Jacobi方法, Gauss-Seidel方法

我们能不能使用类似不动点迭代的思想进行求解呢？我们要考虑：

• 怎么选择迭代形式
• 迭代要收敛
• 收敛的速度也要保证
  

2.1 Jacobi方法

$$\begin{array}{rl}A& =L+D+U\\ x& =D^{-1}(b-Lx-Ux)\\ x^{(k+1)}& =D^{-1}\left(b-Lx{}^{(k)}-{Ux}^{(k)}\right)\end{array}$$

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

例子

具体编程的实现可以有：

或者

然而同一个方程组不同的方程顺序可能会不收敛。

所以什么时候会收敛呢？

定义：
矩阵$A$严格对角占优：$$\left|a_{kk}\right|>\sum _{j=1,j\ne k}^N\left|a_{kj}\right|k=1,2,\dots ,N$$
矩阵$A$表示为：
$$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1N}\\ a_{21}& a_{22}& \cdots & a_{2N}\\ ⋮& ⋮& & ⋮\\ & & & \\ a_{N1}& a_{N2}& \cdots & a_{NN}\end{array}\right]$$

例子：


严格对角占优只能算是Jacobi方法的一个充分条件。

2.2 Gauss-Seidel方法

$$\begin{array}{rl}& A=L+D+U\\ & x=D^{-1}(b-Lx-Ux)\\ & x^{(k+1)}=D^{-1}\left(b-L{x}^{(k+1)}-U{x}^{(k)}\right)\end{array}$$

$$\begin{array}{rl}& A=L+D+U\\ & (L+D)x^{(k+1)}+U{x}^{(k)}=b\\ & x^{(k+1)}=-(L+D{)}^{-1}U{x}^{(k)}+(L+D{)}^{-1}b\end{array}$$

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

迭代速度比Jacobi方法更快！

2.3 Jacobi方法, Gauss-Seidel方法收敛的条件

充分条件1

矩阵$A$满足下面任一①严格对角占优 ②**不可约矩阵（回顾前置知识3），使用Jacobi方法, Gauss-Seidel方法都收敛。

充分必要条件

矩阵的谱半径小于1。
$$\rho (B)=\underset{i}{\max }\left|{\lambda }_i\right|<1\left|B-{\lambda }_{i}I\right|=0$$

以下是对实对称阵说明

$$\begin{array}{rl}& x=Bx+g\\ & x^{(k+1)}=B{x}^{(k)}+g\\ & e^{(k+1)}=B{e}^{(k)}=B^{k+1}{e}^{(0)}【e^{(k+1)}=x^{(k+1)}-x^{(k)}】\\ & e^{(0)}=c_1{V}_1+c_2{V}_2+\dots +c_n{V}_n,B{V}_i={\lambda }_i{V}_i\\ & e^{(k+1)}=c_1{\lambda }_1^{k+1}{V}_1+c_2{\lambda }_2^{k+1}{V}_2+\dots +c_n{\lambda }_n^{k+1}{V}_n\to 0\text{ when }\left|{\lambda }_i\right|<1\end{array}$$

对一般的矩阵形式可以用Jordan形式证明

$$\begin{array}{rl}& B=TJ{T}^{-1}\\ & B^k=T{J}^k{T}^{-1}\\ & J^k=\mathrm{diag}\left(J_{r_1}^k\left({\lambda }_1\right),J_{r_2}^k\left({\lambda }_2\right),\dots ,J_{r_p}^k\left({\lambda }_p\right)\right)\to 0,\text{ when }\left|{\lambda }_i\right|<1\\ & J_{r_i}\left({\lambda }_i\right)=\left[\begin{array}{cccc}{\lambda }_i& 1& & \\ & {\lambda }_i& 1& \\ & & \ddots & 1\\ & & {\lambda }_i& \end{array}\right]\end{array}$$

定理：对于矩阵的任意范数，谱半径都小于矩阵的范数。
$$\rho (B)\le \Vert B\Vert$$

于是有

充分条件2

满足条件$\Vert B\Vert <1$，矩阵$A$使用Jacobi方法, Gauss-Seidel方法都收敛。

说明：以下等价
$$B^k\to 0\iff \Vert B^k\Vert \to 0\iff \rho (B)<1$$

$$\Vert B^k\Vert \le \Vert B\Vert \Vert B^{k-1}\Vert \le \Vert B{\Vert }^k$$

充分条件3
定理：矩阵$A$满足下面任一条件：①严格对角占优 ②**不可约矩阵（回顾前置知识3），并且对角线上元素大于0，$A$是个正定矩阵！

如果矩阵$A$是对称正定矩阵，使用Jacobi方法, Gauss-Seidel方法都收敛。

一张关系图说明收敛的条件

2.4 预测迭代次数

定理：设基本迭代的迭代矩阵$\Vert B\Vert =q<1$,若$\Vert x^{(k+1)}-x^{(k)}\Vert ⩽\epsilon$,则$\Vert x^{(k)}-x\Vert ⩽\frac{\epsilon }{1-q}$

容易证明：$$\begin{array}{r}\Vert X_k-X^{\ast }\Vert \le \frac{1}{1-\Vert B\Vert }\Vert X_{k+1}-X_k\Vert \\ \Vert X_k-X^{\ast }\Vert \le \frac{\Vert B{\Vert }^k}{1-\Vert B\Vert }\Vert X_1-X_0\Vert \end{array}$$

这个定理使用的两种方式：
1.预测需要的迭代次数
2.使用$|x_{k+1}-x_k|$看是否停止迭代。

预测迭代次数类似不动点迭代收敛的推导： 就不具体展开了。

2.5 连续超松弛方法(The Method of Successive Over-Relaxation【SOR】)

$$\begin{array}{rl}& (L+D){\tilde{x}}^{(k+1)}+U{x}^{(k)}=b\\ & x^{(k+1)}=\omega {\tilde{x}}^{(k+1)}+(1-\omega )x^{(k)}\\ & x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\frac{\omega }{a_{ii}}\left(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}\right)\\ & =(1-\omega )x_i^{(k)}+\frac{\omega }{a_{ii}}\left(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-a_{ii}{x}_i^{(k)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}+a_{ii}{x}_i^{(k)}\right)\\ & =x_i^{(k)}+\frac{\omega }{a_{ii}}\left(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i}^n{a}_{ij}{x}_j^{(k)}\right)\\ \Rightarrow & x^{(k+1)}=x^{(k)}+\omega D^{-1}\left(b-L{x}^{(k+1)}-D{x}^{(k)}-U{x}^{(k)}\right)\\ \Rightarrow & D{x}^{(k+1)}=D{x}^{(k)}+\omega \left(b-L{x}^{(k+1)}-D{x}^{(k)}-U{x}^{(k)}\right)\\ \Rightarrow & (D+\omega L)x^{(k+1)}=[(1-\omega )D-\omega U]x^{(k)}+\omega b\\ \Rightarrow & x^{(k+1)}=(D+\omega L{)}^{-1}[(1-\omega )D-\omega U]x^{(k)}+\omega (D+\omega L{)}^{-1}b\end{array}$$

$$x^{(k+1)}=B{x}^{(k)}+c,B=(D+\omega L{)}^{-1}[(1-\omega )D-\omega U],c=\omega (D+\omega L{)}^{-1}b$$

还有另一种形式：
$$\begin{array}{rl}& (L+D){\tilde{x}}^{(k+1)}+U{x}^{(k)}=b\\ & x^{(k+1)}=\omega {\tilde{x}}^{(k+1)}+(1-\omega )x^{(k)}\\ \Rightarrow & x^{(k+1)}=\omega \left(-(L+D{)}^{-1}U{x}^{(k)}+(L+D{)}^{-1}b\right)+(1-\omega )x^{(k)}\\ \Rightarrow & x^{(k+1)}=\left[(1-\omega )I-\omega (L+D{)}^{-1}U\right]x^{(k)}+\omega (L+D{)}^{-1}b\end{array}$$

$$x^{(k+1)}=B{x}^{(k)}+c,B=\left[(1-\omega )I-\omega (L+D{)}^{-1}U\right],c=\omega (L+D{)}^{-1}b$$

$可以根据\omega$取值不同分类如下：

$$\begin{array}{rl}& 0<\omega <2\\ & \omega =1:\text{ Gauss - Seidel }\\ & 0<\omega <1:\text{ Under - Relaxation }\\ & 1<\omega <2:\text{ SOR }\end{array}$$

SOR 收敛的条件和Jacobi方法, Gauss-Seidel方法相同：

1.当系数矩阵A为强对角占优矩阵时，SOR方法收敛；
2.当系数矩阵A为不可约对角占优矩阵时，SOR方法收敛；
3.当系数矩阵A为对称正定矩阵时，SOR方法收敛。

$$x^{(k+1)}=B{x}^{(k)}+c$$

$B$是关于松弛因子$\omega$的一个函数，所以$\omega$应该取多少呢？

可惜的是，${\omega }_{opt}$无法准确求得，只能估算，下面给出两种估算方法。

方法1: 先用$\omega =1$算得$x^{(1)}$和$x^{(2)}$，再用$\omega =1.1$算得${\tilde{x}}^{(1)}$和${\tilde{x}}^{(2)}$;比较$\Vert x^{(1)}-x^{(2)}\Vert$和$\Vert {\tilde{x}}^{(1)}-{\tilde{x}}^{(2)}\Vert$的大小，量值$\Vert {\tilde{x}}^{(1)}-{\tilde{x}}^{(2)}\Vert$较大说明取$\omega =1.1$时迭代收敛快：继续选$\omega =1.2$计算且与$\omega =1.1$的情形比较，不断改进$\omega$的值直到接近$\omega$为止、
方法2: 用$\omega =1.9$和$\omega =1.8$计算，判断比较相应松弛迭代收敛的快慢表现；不断改进参数w的取值，在${\omega }_{opt}$附近还可作些适当的微调处理。

例子：

在SOR方法计算下，当系数矩阵$B$的谱半径小于1但是非常接近1的时候，收敛速度较慢。

来看下面这个例子：

2.6 总结：Jacobi方法, Gauss-Seidel方法对比：

2.7 对称矩阵的Gauss-Seidel方法

对于一个对称矩阵$A$

$$A=L+D+L^T$$

$D$是$A$的对角线组成的对角阵，$L$是$A$的下三角矩阵，由于对称，$L^T$是$A$的上三角矩阵。

令$M=(L+D)D^{-1}(L+D{)}^T$

$$\begin{array}{rl}{x}^{(k+1)}& =x^{(k)}+M^{-1}\left(b-A{x}^{(k)}\right)\\ & =M^{-1}b+M^{-1}(M-A)x^{(k)}\\ & =M^{-1}b+M^{-1}(M-L-D-L^T)x^{(k)}\\ & =M^{-1}b+M^{-1}((L+D)D^{-1}(L+D{)}^T-L-D-L^T)x^{(k)}\\ & =M^{-1}b+M^{-1}(L{D}^{-1}{L}^T+L{D}^{-1}{D}^T+D{D}^{-1}{L}^T+D{D}^{-1}{D}^T-L-D-L^T)x^{(k)}\\ & =M^{-1}b+M^{-1}(L{D}^{-1}{L}^T+L+L^T+D^T-L-D-L^T)x^{(k)}\\ & =M^{-1}b+M^{-1}L{D}^{-1}{L}^T{x}^{(k)}\\ & =M^{-1}b+B{x}^{(k)}\end{array}$$

式中的$B$：

$$\begin{array}{rl}B& =M^{-1}L{D}^{-1}{L}^T\end{array}$$

2.8 Krylov方法(Krylov methods)

$$\begin{array}{rlr}& Ax=b& \\ & q(\lambda )=|A-\lambda I|=a_0+a_1\lambda +\dots +a_n{\lambda }^n& =0\\ & q(A)=a_{0}I+a_{1}A+\dots +a_n{A}^n=0& \\ & -\frac{1}{a_0}A\left(a_{1}I+\dots +a_n{A}^{n-1}\right)=I& \\ & A^{-1}=-\frac{1}{a_0}\left(a_{1}I+\dots +a_n{A}^{n-1}\right)& \\ & x=A^{-1}b=-\frac{1}{a_0}\left(a_{1}b+\dots +a_n{A}^{n-1}b\right)\in \mathrm{span}\left(b,Ab,A^{2}b,\dots ,A^{n-1}b\right)& \\ & x^{\ast }=\sum _i{c}_i{A}^{i}b& \end{array}$$

$x^{\ast }$的维数可小于$n$！

给定

$$\begin{array}{l}{x}_0=0\\ x_n=\left[b,Ab,A^{2}b,\dots ,A^{n-1}b\right]\tilde{c}\end{array}$$

$x_n$和$x^{\ast }$落在一个空间，但可能不是一个很好的近似。

（这个图有点不清楚，后面有空改一下）

原问题的解罗落在于一个Krylov 空间，其维数是$A$的最小多项式的维度。因此，如果 $A$ 的最小多项式的次数较低，则空间维数可以很小。

原来的空间张成向量$b,Ab,A^{2}b,\cdots ,A^{n-1}b$不正交转化为新的标准正交基$q_q,q_2,\cdots ,q_n$就有了下面的方法。

2.9 GMRES方法(Generalized Minimum Residual Method)

解释上面的第二步：

利用QR分解

3. 【线性方程组】基于优化的方法

• 若$A$是一个对称正定矩阵，下面问题等价：

$$Ax=b\iff \underset{x}{\min }\frac{1}{2}{x}^{T}Ax-b^{T}x$$

• 若$A$是一个大型稀疏方阵，下面问题等价：

$$Ax=b\iff \underset{x}{\min }\Vert Ax-b\Vert$$

3.1 共轭梯度法(Conjugate gradient method【CG】)

两个向量$S_1$和$S_2$是共轭的，当它们满足$S_1^{\mathrm{T}}A{S}_2=0$

$$\begin{array}{rl}& f(X)=\frac{1}{2}{X}^{T}AX-b^{T}X\\ & \nabla f(X)=AX-b\\ & {\varphi }^{\prime }\left(a_1\right)=S_i^T\nabla f\left(X_1+a_1{S}_i\right)=S_i^T\left(A{X}^{(1)}-b\right)=0X^{(1)}=X_1+a_1{S}_i\\ & {\varphi }^{\prime }\left(a_2\right)=S_i^T\nabla f\left(X_2+a_2{S}_i\right)=S_i^T\left(A{X}^{(2)}-b\right)=0X^{(2)}=X_2+a_2{S}_i\\ & S_i^{T}A\left(X^{(2)}-X^{(1)}\right)=S_i^{T}AS=0\end{array}$$
可以有限步数收敛！

算法流程：
$$\begin{array}{rl}& 1.g_0=\nabla f\left(x_0\right),d_0=-g_0;\\ & 2.When\left|g_k\right|<eps,exit;\\ & 3.a_k=\underset{a_k}{\mathrm{argmin}}f\left(x_k+a_k{d}_k\right);\\ & 4.g_{k+1}=\nabla f\left(x_{k+1}\right)=\nabla f(x_k+a_k{d}_k);\\ & 5.{\beta }_k=g_{k+1}^T{g}_{k+1}/g_k^T{g}_k;\\ & 6.d_{k+1}=-g_{k+1}+{\beta }_k{d}_k;\\ & 7.k=k+1,goto2.\end{array}$$

共轭梯度法(CG)特点：
①有限次收敛迭代
②计算复杂度$O(n^3)$，如果矩阵$A$是一个对角矩阵，计算复杂度下降为$O(\omega n^2)$

3.2 共轭梯度法和其他算法的对比

共轭梯度法(CG) vs 高斯消元

①CG在所有有限步数之前可以得到足够精确的解
②CG保证稀疏性，即使$A$不是对角阵。

共轭梯度法 vs 迭代方法

①CG保证收敛.
②收敛速度不同。

$${\Vert x^{(k+1)}-x^{\ast }\Vert }_A<C{\Vert x^{(k)}-x^{\ast }\Vert }_A,C=\frac{\sqrt{\mathrm{cond}(A)}-1}{\sqrt{\mathrm{cond}(A)}+1},\Vert x{\Vert }_A=\sqrt{x^{T}Ax}$$

3.3 预条件共轭梯度法(Preconditioned conjugate gradient method)

在1.3我们讲到了条件数，预条件共轭梯度法就是想办法找到合适的$\tilde{A}$降低条件数的大小，提高解的稳定性。

我们这样构造$\tilde{A}$以及$\tilde{x}$和$\tilde{b}$

$$\tilde{A}=C^{-1}A{C}^{-T},\tilde{x}=C^{T}x,\tilde{b}=C^{-1}b$$

于是有：
$$Ax=b\iff \tilde{A}\tilde{x}=\tilde{b}$$

我们可以验证：

$${cond}_2(\tilde{A})<{cond}_2(A){cond}_2(A)=\frac{\left|{\lambda }_{\max }\right|}{\left|{\lambda }_{\min }\right|}$$

1. 按照下面的方式取定$C$，条件数变小了【这一块没写清楚，回头再写写补坑】

$$C^{-1}=D^{-1/2},D=\mathrm{diag}\left(a_{11},a_{22},\dots ,a_{nn}\right),A={\left(a_{ij}\right)}_{n\times n}$$

2. 按照下面的方式取定$C$，条件数变小了

预条件共轭梯度法算法流程：

4. 【非线性方程组】Jacobian矩阵、Newton迭代法、不动点迭代法、Seidel迭代法

4.1 Jacobian矩阵

针对一个非线性方程组，可能有无穷多解。利用Jacobian矩阵迭代线性化处理变为求解线性方程组。

利用Jacobian矩阵引出Newton方法：

4.2 Newton迭代法

算法流程

Jacobian矩阵的问题

①出现奇异，内部元素分母为0或矩阵的秩为0
②收敛阶数降低
③有不确定性

4.3 不动点迭代法

假设 $g_i$ 的偏导数在包含不动点 $(p，q，r)$的一个区域上是连续的。如果选择的起点足够接近定点，并且
$$\begin{array}{l}\left|\frac{\partial g_1}{\partial x}(p,q,r)\right|+\left|\frac{\partial g_1}{\partial y}(p,q,r)\right|+\left|\frac{\partial g_1}{\partial z}(p,q,r)\right|<1,\\ \left|\frac{\partial g_2}{\partial x}(p,q,r)\right|+\left|\frac{\partial g_2}{\partial y}(p,q,r)\right|+\left|\frac{\partial g_2}{\partial z}(p,q,r)\right|<1,\\ \left|\frac{\partial g_3}{\partial x}(p,q,r)\right|+\left|\frac{\partial g_3}{\partial y}(p,q,r)\right|+\left|\frac{\partial g_3}{\partial z}(p,q,r)\right|<1,\end{array}$$

这意味着$\rho (J)\le \Vert J{\Vert }_{\infty }<1$。
那么不动点迭代收敛到定点。

4.4 Seidel迭代法

5. Matlab相关函数