数值分析中的方程组的迭代和直接求解

方程组的迭代&直接求解

一、迭代算法

理论基础-不动点理论

$$x=\varphi (x),$$这样一种映射$\varphi (x)$存在使得$x$的函数值等于其本身，这样的点看起来像是不动点
Lipzchi条件 $$||x_{k+1}-x_k||\le L||\varphi (x_{k+1})-\varphi (x_k)||$$
满足该条件的迭代公式构造的点列收敛

1.一元变量函数的迭代求解

简单迭代

构造迭代函数 $x=\varphi (x)$

各种魔改的牛顿迭代

将$f(x)=0$线性化$L(x)=0$，$L(x)=f(x)+f^{\prime }(x)(x-x\ast )=0$,构造迭代函数$$x^{k+1}=x^k-f(x)/f^{\prime }(x)$$

2.大型方程组（多元变量）的迭代求解

大型方程组构成 $Ax=b$ 的 $m\ast n$ 阶矩阵，若使用直接求解法，例如经典而优美的克莱姆法则，虽然可以得到准确解$X\ast$，但其计算量太高，不适用于大型矩阵的求解；采用迭代法不断逼近真实解$X\ast$

朴素而经典-Jacobi迭代

改进版本的-Seidel迭代

SOR迭代（松弛迭代）

收敛条件

1个充要条件，3个判别条件:
迭代矩阵B的谱半径小于1，
$$\rho (B)<1$$
$det(\lambda I-B)=0$ 解的最大值 $\lambda max$

范数

衡量某种距离$$||\cdot ||$$， 向量范数与矩阵范数不同
常见向量范数有： 1范数，2范数，p范数，无穷范数

常见矩阵范数有： 1范数（列范数），无穷范数（行范数），F范数（元素平方和），2范数（乘转置矩阵求特征值）

现代大规模矩阵迭代求解方法

参考资料 Yosse Sparse Matrix iterative method
比较常用的是Krylov空间方法
GMRes （Generalization minimize residual），CG （conjugate gradient），…

先看看对 $Ax=b$ 的理解

1. 朴素线性方程组
   $a_{11}{X}_1+a_{12}{X}_2+a_{13}{X}_3+...+a_{1n}{X}_n=0$
   $a_{21}{X}_1+a_{22}{X}_2+a_{23}{X}_3+...+a_{2n}{X}_n=0$
   …
   这里的系数矩阵 $A$=
   { $a_{11}$ $a_{12}$ $a_{13}$ …$a_{1n}$
   $a_{21}$ …
   … }
   可以从一个微分方程组，空间离散/线性化后的得到的；

2. 牛顿Newton角度
   让我们看看牛顿的问题
   求解非线性方程 $x^n+x^{n-1}+x^{n-2}+.....+x^1=b$
   使用 $f(x)=x^n+x^{n-1}+x^{n-2}+.....+x^1-b$
   则问题转化为 $f(x)=0$ 的求根问题
   典型的牛顿Newton迭代解法
   $x^{k+1}=x^k-f(x^k)/f^{\prime }(x^k)$，$f^{\prime }(x^k)$是原函数$f(x)$在$x^k$的导数，对于多项式非线性这很好计算；
   ok，对于多个非线性方程组则可以构成一个$Ax=b$的形式；
   $x_1^n+x_2^{n-1}+x_3^{n-2}+.....+x_n^1=b_1$
   $ex{p}^{x_1}+x_2+x_3^n+.....+x_n^1=b_2$
   …
   其变量 $\mathbf{x}=(x_1,x_2,x_3,...{)}^T$ 是一个向量vector，自然其函数 $\mathbf{f}=(f_1,f_2,f_3,...)$ 也构成一个向量vector。
   则其迭代计算也变成了
   ${\mathbf{x}}^{\mathbf{k}+1}={\mathbf{x}}^{\mathbf{k}}-\mathbf{f}({\mathbf{x}}^{\mathbf{k}})/{\mathbf{f}}^{\prime }({\mathbf{x}}^{\mathbf{k}})$
   这里${\mathbf{f}}^{\prime }({\mathbf{x}}^{\mathbf{k}})$构成一个Jacob矩阵 $\mathbf{J}$=
   $\partial f_1/\partial x_1,\partial f_1/\partial x_2,....$
   $\partial f_2/\partial x_1,\partial f_2/\partial x_2,....$
   …
   则 ${\mathbf{x}}^{\mathbf{k}+1}={\mathbf{x}}^{\mathbf{k}}-{\mathbf{J}}^{-1}\mathbf{f}({\mathbf{x}}^{\mathbf{k}})$
   令 $\delta {\mathbf{x}}^k={\mathbf{x}}^{k+1}-{\mathbf{x}}^k$，再左乘$\mathbf{J}$得到
   $\mathbf{J}$$\delta {\mathbf{x}}^k$=$\mathbf{f}({\mathbf{x}}^{\mathbf{k}})$
   即 $Ax=b$形式，其中$A$为Jacobain矩阵$\mathbf{J}$, $b$右端向量为$\mathbf{f}({\mathbf{x}}^{\mathbf{k}})$。
   这二者每步迭代都需要在${\mathbf{x}}^{\mathbf{k}}$处更新。

3. 最优化凸优化角度
   定义 $r=Ax-b$，这里$r$是余向量 redisual vector，在 $r=0$ 的时候$Ax=b$，所以迭代问题转化为最小化这个余向量 $r$，这个$r$的本质是最优化问题$f(x)=x^{T}Ax-x^{T}b$的梯度方向向量。
   可以简单证明，若$A$为对称正定矩阵（SPD），则 $f(x)=x^{T}Ax-x^{T}b$ 的导数为$Ax-b$,导数为零取得优化解，即满足 $r=Ax-b=0$
   GMRes中，迭代按照搜索方向$p^k$进行,$p^k$与余向量 $r$和A都有关系，力求按照$\alpha$步长每次在该方向使得$||r||$到达极小 即minizie。
   例如Krylov家族的CG方法：迭代公式为 $x^{k+1}=x^k+{\alpha }^k{p}^k$
   此方法最多需要n步即可得出解，远快于朴素的迭代方法$x=\varphi (x)$，很适合大规模矩阵求解。

   GMRes包括Arnoldi求正交基步骤（对Krylov空间求正交基）；Gevins 旋转步骤
   Krylov空间即向量 $(b,Ab,A^{2}b,...A^{n-1}b)$span 张成的空间，其向量构成的矩阵即Krylov矩阵，各向量不一定正交，不一定正定。

3. 方程组的直接求解

Gauss消元法

实质是将$Ax=b$中的A矩阵转化为同等的上三角矩阵$U$,然后回代，解出$x$

LU分解

将A矩阵分解为两个上下三角矩阵LU的乘积
$$A=LU$$ $$Ax=b$$ $$LUx=b$$ $$Ly=b,Ux=y$$

Doolittle分解
紧凑格式

LDU分解

$$A=LDU$$

追赶法

解决三对角矩阵，即$A$是稀疏的，A呈现条带状分布
计算量少： 5n-4 + 3n-2 乘除法
存储简单：一维数组存储，a，b，c以及 alpha, beta, gamma

关注问题：大规模稀疏方程组

（1）三对角矩阵可以用直接求解：追赶法，Guass-Doolittle

（2）分块迭代

（3）CG共轭梯度法,共轭方向法的一种，利用当前梯度方向构建共轭方向
所谓共轭向量组是满足 $$d_{i}A{d}_j=0$$的向量组,其中A是正定矩阵，$x^{T}Ax$>0
$S=d_0，d_1,d_2,...d_k$,共轭向量组$S$是线性无关的，而且是正交的
其沿着共轭方向$d_i$,通过$$x_{k+1}=x_k+{\alpha }_k{d}_k$$构造点列$x_0,x_1,x_2,...,x_n$,
优势：具有n步平方收敛
理论基础是 最优化中的求极值？而不是不动点理论

关注问题：矩阵病态问题

误差被放大倍数与条件数 cond(A)相关 $$cond(A)=||A||\cdot ||A^{-1}||$$
cond(A)>>1,矩阵是病态的
矩阵预处理，降低条件数cond(A)
利用剩余向量 $$R=Ax-b，$$精确结果,直到$||R||<ϵ$