GAMES001课程笔记12-线性系统

这个系列是GAMES001-图形学中的数学(GAMES 001: Mathematics in Computer Graphics)的同步课程笔记。课程旨在总结归纳图形学学习过程中重要的数学概念、理论和方法，并将以索引的形式在每一章节中穿插讲解该数学专题在图形学中的应用。本课程既可以作为GAMES系列其它课程学习的基础或「手册」，也可以作为站在不一样的视角复习图形学知识的平台。本节主要介绍求解线性系统的相关技术。

求解线性系统实际上就是解矩阵方程

\[\mathbf{A} \boldsymbol{x} = \mathbf{b}\]

在图形学的各个研究方向中，许多问题的本质就是求解一个线性系统。

直接求解

高斯消元

求解线性系统的基本方法是直接进行求解，其中高斯消元法是这类方法的典型代表。具体来说，高斯消元法会同时对系数矩阵\(\mathbb{A}\)和向量\(\mathbf{b}\)进行行变换，使得矩阵\(\mathbb{A}\)变成一个上三角矩阵。得到上三角矩阵后即可从下往上反向回代计算出解\(\boldsymbol{x}\)。

换个角度来看，矩阵的行变换等价于左乘一个下三角矩阵。因此，高斯消元法可以理解为对矩阵进行LU分解。通过LU分解，我们将系数矩阵\(\mathbf{A}\)分解为一个下三角矩阵\(\mathbf{L}\)和一个上三角矩阵\(\mathbf{U}\)的乘积。然后，通过求解两个较简单的线性系统\(\mathbf{L} \mathbf{y} = \mathbf{b}\)和\(\mathbf{U} \boldsymbol{x} = \mathbf{y}\)就可以得到最终解\(\boldsymbol{x}\)。

从编程的角度来讲，高斯消元法的核心在于计算矩阵的LU分解，其伪代码可以参考下方ppt。显然LU分解的复杂度为\(O(n^3)\)，因此LU分解是一个复杂度比较高的算法，不过可以利用一些并行的方法来进行加速。

使用高斯消元法时还需要注意，如果当前选择的对角线元素值接近0时会产生数值计算上的不稳定。为了处理这种情况，可以考虑每次都选择当前列中绝对值最大的那一行作为主元，然后再利用主元来消去其它行。这种带选择主元的高斯消元法也称为列主元高斯消元法或部分主元高斯消元法(Gaussian Elimination with Partial Pivoting, GEPP)。

Cholesky分解

对于对称半正定的系数矩阵，LU分解可以简化为Cholesky分解。此时，系数矩阵\(\mathbf{A}\)可以表示为下三角矩阵\(\mathbf{A}\)与其自身转置\(\mathbf{A}^T\)的乘积。

Cholesky分解的伪代码可以参考下方ppt。类似于LU分解，Cholesky分解的复杂度同样是\(O(n^3)\)，也可以利用一些并行的技术来进行加速。

简单总结一下，LU分解以及高斯消元法可以求解任意的矩阵方程，而Cholesky分解则只适用于对称半正定的系统。由于这两种方法都具有\(O(n^3)\)的复杂度，它们只适合一些规模比较小的线性系统。

迭代求解

对于更大规模的线性系统，我们通常会使用迭代法来求解。迭代法的基本框架是从一个初始解\(x^0\)出发，反复执行某个步骤来更新当前解

\[x^{k+1} \leftarrow \Psi (x^k)\]

迭代过程会继续直到满足设定的迭代次数上限，或者达到一定的残差收敛条件。相较于直接法，迭代法的优势在于可以在任意迭代步骤中得到精度满足要求的近似解，而不需要等到整个算法结束才能终止。在许多应用场景中，仅需一个近似解即可满足需求。同时，一个好的初始解\(x^0\)往往可以极大地减少迭代次数。因此对于特定问题优化的迭代算法有着远高于直接求解的效率。除此之外，迭代法甚至可以推广到系数矩阵\(\mathbf{A}\)没有显式给出的情况。

不动点迭代

首先我们来介绍不动点迭代。举个简单的例子，假设要求解方程

\[x = \cos{x}\]

我们可以建立迭代格式

\[x^{k+1} \leftarrow \cos{x^k}\]

通过不断执行可以发现\(\boldsymbol{x}\)会收敛到某个值上，而实际上这个值就是方程的解。

在线性系统的求解中，可以利用不动点迭代方法来进行计算。我们可以将系数矩阵\(\mathbf{A}\)表示为

\[\mathbf{A} = \mathbf{M} - \mathbf{N}\]

这样线性系统的解可以表示为

\[\boldsymbol{x} = \mathbf{M}^{-1} (\mathbf{N} \boldsymbol{x} + \mathbf{b})\]

从而建立迭代格式

\[\boldsymbol{x}^{k+1} \leftarrow \mathbf{M}^{-1} (\mathbf{N} \boldsymbol{x}^k + \mathbf{b})\]

如果迭代能够收敛，则可以保证\(\boldsymbol{x}^k\)会趋向于线性系统的解。以上就是不动点迭代求解线性系统的基本思想，当然使用不动点迭代时还有很多问题需要考虑，比如如何计算逆阵\(\mathbf{M}^{-1}\)、如何保证迭代收敛以及迭代的效率等。

Jacobi迭代

Jacobi迭代是一种典型的不动点迭代法。在Jacobi迭代中，我们取\(\mathbf{M}\)矩阵为系数矩阵\(\mathbf{A}\)的对角元素\(\mathbf{D}\)，则\(\mathbf{N}\)矩阵为系数矩阵\(\mathbf{A}\)的其余非对角元素的相反数。此时\(\mathbf{M}^{-1}\)为对角元素的倒数，因此每次迭代都只需要进行基本的矩阵乘法和加法运算。

Jacobi迭代的伪代码可以参考如下。

需要注意的是，在某些情况下Jacobi迭代可能会出现不收敛的情况。

要分析不动点迭代法的收敛准则，我们需要引入「误差」的概念。这里「误差」是指当前解\(\boldsymbol{x}^k\)与真实解\(\boldsymbol{x}^*\)之间的差。误差\(\boldsymbol{e}^k\)与残差\(\boldsymbol{r}^k\)之间的关系为

\[\boldsymbol{r}^k = \mathbf{b} - \mathbf{A} \boldsymbol{x}^k = \mathbf{A} (\boldsymbol{x}^* - \boldsymbol{x}^k) = \mathbf{A} \boldsymbol{e}^k\]

显然，随着迭代的进行误差\(\boldsymbol{e}^k\)具有递推关系

\[\boldsymbol{e}^{k+1} \leftarrow (\mathbf{I} - \mathbf{M}^{-1} \mathbf{A}) \boldsymbol{e}^k = \mathbf{M}^{-1} \mathbf{N} \boldsymbol{e}^k\]

上式说明，在每一次迭代时误差更新相当于左乘一个矩阵\(\mathbf{T}\)

\[\mathbf{T} = \mathbf{I} - \mathbf{M}^{-1} \mathbf{A} = \mathbf{M}^{-1} \mathbf{N}\]

只有当矩阵\(\mathbf{T}\)的谱半径小于1时误差才能保证收敛，否则误差可能会不收敛。

在Jacobi迭代中，如果系数矩阵\(\mathbf{A}\)是严格对角占优的则可以保证迭代一定可以收敛。

如果系数矩阵是非对角占优的，还可以通过松弛的方法使得Jacobi迭代能够收敛。

总结来说，Jacobi迭代是一种简单实现且适合并行计算的迭代求解算法。然而，实践经验表明Jacobi迭代通常需要较多的迭代步数才能达到收敛。当系数矩阵接近对角占优时，Jacobi迭代的收敛速度会明显加快。

Gauss-Seidel迭代

除了Jacobi迭代，另一种常用的迭代算法是Gauss-Seidel迭代。和Jacobi迭代相比，Gauss-Seidel迭代在构造\(\mathbf{M}\)矩阵时使用了系数矩阵的对角元以及下三角部分，此时计算\(\mathbf{M}^{-1}\)可以使用类似与高斯消元法的方式进行处理。由于这样构造的\(\mathbf{M}\)矩阵更接近与系数矩阵\(\mathbf{A}\)，因此Gauss-Seidel迭代通常具有比Jacobi迭代更快的收敛速度和更高的收敛效率。特别值得注意的是，Gauss-Seidel迭代的收敛充分但不是必要条件是系数矩阵\(\mathbf{A}\)为对称半正定矩阵。

从代码实现的角度来看，Gauss-Seidel迭代和Jacobi迭代的唯一区别在于更新当前解\(\boldsymbol{x}^{k+1}\)时是否使用已经更新过的结果。Gauss-Seidel迭代会利用当前步中已经更新过的计算结果，而Jacobi迭代只会利用上一步的结果。

两种不动点迭代方法的优劣可以总结如下。

子空间迭代

除了前面介绍的不动点迭代方法，另一种常用的迭代求解线性系统的方法为子空间迭代法。

子空间迭代法中最常用的方法是Krylov子空间迭代法。

在Krylov子空间的基础上，对于对称半正定系数矩阵的情况我们可以推导共轭梯度下降法。其基本思路是将求解矩阵方程\(\mathbf{A} \boldsymbol{x} = \mathbf{b}\)转化为一个最优化问题：

\[\min f(\boldsymbol{x}) = \frac{1}{2} \boldsymbol{x}^T \mathbf{A} \boldsymbol{x} - \mathbf{b}^T \boldsymbol{x}\]

而优化目标函数的梯度为

\[\nabla f = \mathbf{A} \boldsymbol{x} - \mathbf{b} = - \boldsymbol{r}\]

共轭梯度法的核心思想是在每次迭代中，将当前的搜索方向与之前的搜索方向正交化(即共轭)，从而在Krylov子空间 \(K_m\)中进行优化。这个子空间是由初始残差和迭代过程中生成的搜索方向构成的。假设我们能够构造出Krylov子空间的一组基\(K_m = \text{span} \{ \boldsymbol{p}_1, ..., \boldsymbol{p}_m \}\)，它们满足共轭条件：

\[\boldsymbol{p}_i^T \mathbf{A} \boldsymbol{p}_j = 0\]

则在每一步更新\(\boldsymbol{x}_m\)时只需在当前的基上进行移动：

\[\boldsymbol{x}_m = \boldsymbol{x}_{m-1} + \alpha_m \boldsymbol{p}_m\]

其中步长\(\alpha_m\)可以显式求解。可以证明，只要能够构造出一组共轭基\(K_m = \text{span} \{ \boldsymbol{p}_1, ..., \boldsymbol{p}_m \}\)，共轭梯度法是可以运行下去的。