计算方法03 线性方程组求解
Published:
\newcommand\norm[1]{\left\lVert#1\right\rVert} \newcommand\abs[1]{\left\lvert#1\right\rvert}
前置
内积
二元函数 $\left<,\right>$ 被称为内积，则其满足：
$\left<x,y\right>=\left<y,x\right>$
$\left<ax+by,z\right>=a\left<x,z\right>+b\left<y,z\right>$
$\left<x,x\right>\geq 0$，等号仅在 $x=\textbf{0}$ 时取到
向量范数
一元函数 $\norm{\cdot}$ 被称为范数，则其满足：
$\norm{x}\geq 0$，等号仅在 $x=\textbf0$ 时取到
$\norm{kx}=\abs{k}\norm{x}$
$\norm{x+y}\leq \norm{x}+\norm{y}$
在内积空间上有一个天然的范数 $\norm{x}=\sqrt{\left<x,x\right>}$，容易验证满足上述三条要求。
矩阵范数
算子范数
矩阵可以看成线性变换，因此可以由向量范数来衡量矩阵作为线性变换（算子）的“长度”。定义为
\begin{aligned}
\norm{A}=\max_{\norm{x}=1} \norm{Ax}
\end{aligned}
注意这里的 $Ax,x$ 都是向量，因此 RHS 出现的全都是向量范数，LHS 则是定义出来的算子范数。
有了这个定义，我们就可以写出这样的不等式
\begin{aligned}
\norm{Ax}\leq \norm{A}\norm{x}
\end{aligned}
同时会引入新的定义，称满足下述要求的算子范数具有相容性：
\begin{aligned}
\norm{AB}\le \norm{A}\norm{B}
\end{aligned}
有了相容性同样可以做一些界的估计
矩阵范数
矩阵同样也可以看成线性空间中的元素，因此可以单独赋予范数的定义，当然这就和向量范数没什么关系了。
一个例子是这样的，容易验证其满足三条要求
\begin{aligned}
\norm{A}=\sum_{i,j}\abs{A_{i,j} }
\end{aligned}
高斯消元
对于给定的线性方程组 $Ax=b$，可以用简单 $O(n^3)$ 的高斯消元法来求解。并且这样可以很容易地求出解空间的一组基，进而得到所有解。
解的稳定性
不妨假设 $A=xb$ 中 $\abs{A}\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{\norm{A^{-1}\hat b} }{\norm{A^{-1}b} } }/{\frac{\norm{\hat b} }{\norm{b} } }
\end{aligned}
即后向相对误差与前向相对误差的比值，称这个值为方程组 $A$（矩阵 $A$）的条件数 $\text{cond}(A)$
\begin{aligned}
\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}
\end{aligned}
注意到 $A^{-1}b=x$，且 $b=Ax$，带入即得
\begin{aligned}
\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}
\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$ 的谱半径为其绝对值最大的特征值 $\abs{\lambda_\max}$，记为 $\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$ 是个右逆….