计算方法3 线性方程组求解

前置

这一节主要是玩矩阵，为了偷懒只讨论实线性空间

内积

二元函数 \(\left<,\right>\) 被称为内积，则其满足：

1. \(\left<x,y\right>=\left<y,x\right>\)
2. \(\left<ax+by,z\right>=a\left<x,z\right>+b\left<y,z\right>\)
3. \(\left<x,x\right>\geq 0\)，等号仅在 \(x=\bold{0}\) 时取到

向量范数

一元函数 \(\norm{\cdot}\) 被称为范数，则其满足：

1. \(\norm{x}\geq 0\)，等号仅在 \(x=\bold0\) 时取到
2. \(\norm{kx}=\abs{k}\norm{x}\)
3. \(\norm{x+y}\leq \norm{x}+\norm{y}\)

在内积空间上有一个天然的范数 \(\norm{x}=\sqrt{\left<x,x\right>}\)，容易验证满足上述三条要求。

矩阵范数

算子范数

矩阵可以看成线性变换，因此可以由向量范数来衡量矩阵作为线性变换（算子）的“长度”。定义为

\[\norm{A}=\max_{\norm{x}=1} \norm{Ax} \]

注意这里的 \(Ax,x\) 都是向量，因此 RHS 出现的全都是向量范数，LHS 则是定义出来的算子范数。

有了这个定义，我们就可以写出这样的不等式

\[\norm{Ax}\leq \norm{A}\norm{x} \]

同时会引入新的定义，称满足下述要求的算子范数具有相容性：

\[\norm{AB}\le \norm{A}\norm{B} \]

有了相容性同样可以做一些界的估计

矩阵范数

矩阵同样也可以看成线性空间中的元素，因此可以单独赋予范数的定义，当然这就和向量范数没什么关系了。

一个例子是这样的，容易验证其满足三条要求

\[\norm{A}=\sum_{i,j}\abs{A_{i,j}} \]

高斯消元

对于给定的线性方程组 \(Ax=b\)，可以用简单 \(O(n^3)\) 的高斯消元法来求解。并且这样可以很容易地求出解空间的一组基，进而得到所有解。

解的稳定性

不妨假设 \(A=xb\) 中 \(\abs{A}\neq0\)，则有 \(x=A^{-1}b\)

通常 \(A\) 是给定的，而 \(b\) 是若干计算和观察的结果，因此解的误差主要来源于 \(b\) 引入的误差，可以写成

\[\hat x=A^{-1}\hat b \]

考虑相对误差的计算，即为

\[\max_{\hat b,b\neq0} {\frac{\norm{A^{-1}\hat b}}{\norm{A^{-1}b}}}/{\frac{\norm{\hat b}}{\norm{b}}} \]

即后向相对误差与前向相对误差的比值，称这个值为方程组 \(A\)（矩阵 \(A\)）的条件数 \(\text{cond}(A)\)

\[\begin{aligned} \text{cond}(A)&=\max_{\hat b,b\neq 0} {\frac{\norm{A^{-1}\hat b}}{\norm{A^{-1}b}}}/{\frac{\norm{\hat b}}{\norm{b}}} \\ &=\max_{\hat b\neq 0} \frac{\norm{A^{-1}\hat b}}{\norm{\hat b}}\max_{b\neq 0}\frac{\norm{b}}{\norm{A^{-1}b}} \end{aligned} \]

注意到 \(A^{-1}b=x\)，且 \(b=Ax\)，带入即得

\[\begin{aligned} \text{cond}(A)&=\max_{\hat b\neq 0}\frac{\norm{A^{-1}\hat b}}{\norm{\hat b}}\max_{x\neq 0} \frac{\norm{Ax}}{\norm{x}} \\ &=\norm{A^{-1}}\norm{A} \end{aligned} \]

矩阵的条件数反映了矩阵所对应线性方程组的不稳定程度。条件数越大则不稳定程度越大，误差的传递放大也就越严重。

迭代算法

\(O(n^3)\) 太昂贵，考虑迭代算法。一般而言迭代的次数是人为规定的，不少算法能够保证在 \(\Omega(n)\) 次迭代之后必然得到精确解。

Jacobi 迭代

对于矩阵 \(A\)，将其分解为 \(A=L+D+U\)，其中 \(L,U\) 分别为上三角矩阵和下三角矩阵，\(D\) 为对角阵。

那么有

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

收敛性

给出一个充分条件：若 \(A\) 是主对角线占优矩阵，则迭代必然收敛。

只需联立如下方程

\[\left\{ \begin{aligned} x^*&=D^{-1}(b-(L+U)x^*) \\ x_{k+1}&=D^{-1}(b-(L+U)x_{k}) \end{aligned} \right. \]

两式相减即得

\[x_{k+1}-x^*=-D^{-1}(L+U)(x_k-x^*) \]

记 \(W=D^{-1}(L+U)\)，由对角占优可知 \(Wx\) 的每一项绝对值严格小于 \(x\)，因此迭代收敛。

正确性

假设收敛，则有不动点 \(x\)，简单替换即得不动点 \(x\) 是原方程的一个解。

结合收敛性的充分条件即得：若 \(A\) 为主对角线占优矩阵，则解唯一（\(A\) 可逆），且迭代收敛至唯一解。

看起来还是比较好的。

Gauss-Seidel 迭代

用到了一个观察，即在计算解向量 \(x_{k+1}\) 的第 \(r\) 项时，它的前 \(r-1\) 项都已经算出来了（废话）

因此可以写成

\[x_{k+1}=D^{-1}(b-Ux_k-Lx_{k+1}) \]

证明和上面是类似的，结论也是类似的。

写代码可以发现这两个方法没有绝对的好坏之分，迭代速度也没有一般性的结论（至少俺不知道）。

谱半径

为了分析一类迭代算法的收敛性，引入谱半径的概念

定义 \(A\) 的谱半径为其绝对值最大的特征值 \(\abs{\lambda_\max}\)，记为 \(\rho(A)\)

对于这样的迭代算法

\[x_{k+1}=Ax_k+b \]

取不动点做差得

\[x_k-x^*=A^k(x_0-x^*) \]

如果能说明 \(\lim\limits_{k\to\infty} A^k=O\)，那么就能说明迭代是收敛的，且收敛到方程组的解。

有定理：\(\lim\limits_{k\to\infty}A^k=O\) 当且仅当 \(\rho(A)<1\)

分析 Jacobi 和 Gauss-Seidel 迭代矩阵的谱半径可以知道，它们以任意初始向量开始迭代都会收敛。