计算方法03 线性方程组求解
2022-05-08 14:33:00

前置

内积

二元函数 \(\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=\textbf{0}\) 时取到

向量范数

一元函数 \(\left\lVert\cdot\right\rVert\) 被称为范数,则其满足:

  1. \(\left\lVert x\right\rVert\geq 0\),等号仅在 \(x=\textbf0\) 时取到

  2. \(\left\lVert kx\right\rVert=\left\lvert k\right\rvert\left\lVert x\right\rVert\)

  3. \(\left\lVert x+y\right\rVert\leq \left\lVert x\right\rVert+\left\lVert y\right\rVert\)

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

矩阵范数

算子范数

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

\[\begin{aligned} \left\lVert A\right\rVert=\max_{\left\lVert x\right\rVert=1} \left\lVert Ax\right\rVert \end{aligned}\]

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

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

\[\begin{aligned} \left\lVert Ax\right\rVert\leq \left\lVert A\right\rVert\left\lVert x\right\rVert \end{aligned}\]

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

\[\begin{aligned} \left\lVert AB\right\rVert\le \left\lVert A\right\rVert\left\lVert B\right\rVert \end{aligned}\]

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

矩阵范数

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

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

\[\begin{aligned} \left\lVert A\right\rVert=\sum_{i,j}\left\lvert A_{i,j} \right\rvert \end{aligned}\]

高斯消元

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

解的稳定性

不妨假设 \(A=xb\)\(\left\lvert A\right\rvert\neq0\),则有 \(x=A^{-1}b\)

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

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

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

\[\begin{aligned} \max_{\hat b,b\neq0} {\frac{\left\lVert A^{-1}\hat b\right\rVert }{\left\lVert A^{-1}b\right\rVert } }/{\frac{\left\lVert\hat b\right\rVert }{\left\lVert b\right\rVert } } \end{aligned}\]

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

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

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

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

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

迭代算法

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

Jacobi 迭代

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

那么有

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

收敛性

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

只需联立如下方程

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

两式相减即得

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

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

正确性

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

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

看起来还是比较好的。

Gauss-Seidel 迭代

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

因此可以写成

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

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

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

谱半径

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

定义 \(A\) 的谱半径为其绝对值最大的特征值 \(\left\lvert\lambda_\max\right\rvert\),记为 \(\rho(A)\)

对于这样的迭代算法

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

取不动点做差得

\[\begin{aligned} x_k-x^*=A^k(x_0-x^*) \end{aligned}\]

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

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

分析 Jacobi 和 Gauss-Seidel 迭代矩阵的谱半径可以知道,当系数矩阵 \(A\) 是严格对角占优时,这两个算法以任意初始向量开始迭代都会收敛。

伪逆

对于满秩矩阵 \(A\),方程组 \(Ax=b\) 的解是显然存在的。但是对于非满秩的矩阵 \(B\) 来说就不一定了。在最小二乘法中可以求得一个二范数最小的“逼近”解,实际上是解了一个方程 \(A^\intercal Ax=A^\intercal b\),即 \(x=(A^\intercal A)^{-1}A^\intercal b\)

问题的关键在于 \(A^{-1}\) 不一定存在,这使得 \(Ax=b\) 不一定存在唯一解。因此引入记号 \(A^\dagger=(A^\intercal A)^{-1}A^\intercal\),称为 \(A\) 的伪逆(pseudo inverse),那么最小二乘法可以写成 \(x=A^\dagger b\),形式与满秩方程组 \(x=A^{-1}b\) 是一致的。这样求出的解为最佳平方逼近解。

\(A\) 的列线性无关时,\(A^\intercal A\) 满秩,\(A^\dagger=(A^\intercal A)^{-1}A^\intercal\) 存在。此时原方程组 \(Ax=b\) 无解(超定方程组)

由反证法假设存在 \(x\neq 0\) 使得 \(A^\intercal Ax=0\),那么 \((Ax)^\intercal (Ax)=0\),根据内积的正定性可知 \(Ax=0\),这说明存在列的线性组合为 \(0\),这与列线性无关矛盾;故假设不成立。

类似取 \(A^\intercal\) 可知,当 \(A\) 的行线性无关时,\(A^\intercal\) 列线性无关,\(AA^\intercal\) 满秩,\(A^\dagger=((A^\intercal)^\dagger)^\intercal=A^\intercal(AA^\intercal)^{-1}\)。此时方程组 \(Ax=b\) 有多解(欠定方程组),但是我不知道有啥意义,因为 \(A^\dagger\) 是个右逆....