【GAMES103学习笔记】线性系统（Linear System）及其解法（Linear Sovler）

0. 为什么要研究线性系统
1. 什么是线性系统
2. 解线性系统
- 2.1 直接法
- - 2.1.1 直接法的优劣
- 2.2 迭代法
- - 2.2.1 收敛性证明
  - 2.2.2 关于迭代矩阵
  - 2.2.3 迭代法的优劣

小破站链接：https://www.bilibili.com/video/BV12Q4y1S73g?p=3
课程主页：http://games-cn.org/games103/

0. 为什么要研究线性系统

计算机图形学中，很多计算问题，变换到最终都会变换到一个线性系统，以求得最终的结果。

1. 什么是线性系统

对线性系统最直观的理解，就是我们小时候学习的多元一次线性方程组。

10 x + 3 y + 7 z = 2 3 x + 13 y + 12 z = 4 \begin{alignedat}{3} 10x&+& 3y&+& 7&z = 2 \\ 3x&+&13y&+&12&z = 4 \end{alignedat} 10x3x++3y13y++712z=2z=4

但是学习了线性代数后，我们要用矩阵+向量的数学思维和符号去描述、理解线性系统。

一个线性系统通常写作以下形式：

A x = b \bf Ax=b Ax=b

其中， A \bold A A是一个m行n列的矩阵， x \bold x x是一个n元向量， b \bold b b是m元向量。

如上例的线性方程组改写成线性代数形式为：

A = [ 10 3 7 3 13 12 ] , x = [ x y z ] , b = [ 2 4 ] \bold A = \begin{bmatrix} 10 & 3 & 7 \\ 3 & 13 &12 \end{bmatrix}, \bold x = \begin{bmatrix} x\\ y\\ z \end{bmatrix}, \bold b = \begin{bmatrix} 2\\ 4 \end{bmatrix} A=[103313712],x=⎣ ⎡xyz⎦ ⎤,b=[24]

对于矩阵 A \bold A A， m = 2 , n = 3 m=2, n = 3 m=2,n=3， m m m正好是方程的个数， n n n是未知数的个数。

在线形系统中， b \bold b b通常叫做叫做边界向量（boundary vector）。

2. 解线性系统

虽然可以通过求解 A \bold A A的逆，直接求解线性系统。

但是有些情况下 A \bold A A的逆可能不存在，或者 A \bold A A是又巨大、又稀疏（即有好多的0）的方阵，求解耗时、耗空间（对于计算机来说），例如10万个点，每个点都有 x y z xyz xyz坐标值。

我们要寻求更优的解法。

通常有两个常用的方法：直接法（Direct Method）和迭代法（Iterative Method）。

以下内容均假设 A \bold A A是方阵，即 m = n m = n m=n。

2.1 直接法

直接法一般基于矩阵的LU分解（ L U \bf LU LU Factorization），也有很多变种，例如Cholesky、 L D L T LDL^T LDLT等。

A = L U = [ l 00 l 10 l 11 l 20 l 21 l 22 ∣ ⋯ ⋯ ∖ ] [ ∖ ⋯ ⋯ ∣ u n − 2 , n − 2 u n − 2 , n − 1 u n − 2 , n u n − 1 , n − 1 u n − 1 , n u n , n ] \bold A = \bold L \bold U = \begin{bmatrix} l_{00} & \\ l_{10} & l_{11} \\ l_{20} & l_{21} & l_{22}\\ \lvert & \cdots & \cdots & \smallsetminus \end{bmatrix} \begin{bmatrix} \smallsetminus & \cdots & \cdots& \lvert\\ & u_{n-2,n-2} & u_{n-2,n-1} & u_{n-2,n} \\ & & u_{n-1,n-1} & u_{n-1,n} \\ & & & u_{n,n} \end{bmatrix} A=LU=⎣ ⎡l00l10l20∣l11l21⋯l22⋯∖⎦ ⎤⎣ ⎡∖⋯un−2,n−2⋯un−2,n−1un−1,n−1∣un−2,nun−1,nun,n⎦ ⎤

L U \bf LU LU分别是下三角和上三角矩阵。如何得到 L U \bf LU LU本文略过。

得到了 L \bold L L和 U \bold U U之后，首先解线性系统 L y = b \bf Ly=b Ly=b。

因为 L \bold L L是一个下三角矩阵，所以

[ l 00 l 10 l 11 l 20 l 21 l 22 ∣ ⋯ ⋯ ∖ ] [ y 0 y 1 ∣ y n ] = [ b 0 b 1 ∣ b n ] y 0 = b 0 / l 00 y 1 = ( b 1 − l 10 y 0 ) / l 11 ⋯ \begin{bmatrix} l_{00} & \\ l_{10} & l_{11} \\ l_{20} & l_{21} & l_{22}\\ \lvert & \cdots & \cdots & \smallsetminus \end{bmatrix} \begin{bmatrix} y_{0}\\ y_{1}\\ \lvert\\ y_{n} \end{bmatrix} = \begin{bmatrix} b_{0}\\ b_{1}\\ \lvert\\ b_{n} \end{bmatrix}\\ \space \\ \space \begin{aligned} y_{0} &= b_{0} / l_{00} \\ y_{1} &= (b_{1}-l_{10}y_{0}) / l_{11}\\ & \cdots \end{aligned} ⎣ ⎡l00l10l20∣l11l21⋯l22⋯∖⎦ ⎤⎣ ⎡y0y1∣yn⎦ ⎤=⎣ ⎡b0b1∣bn⎦ ⎤ y0y1=b0/l00=(b1−l10y0)/l11⋯

这样一行一行遍历即可把 y \bold y y的值求出来。

之后解线性系统 U x = y \bf Ux=y Ux=y。

因为 U \bold U U是一个上三角矩阵，我们从下往上回解，可得

[ ∖ ⋯ ⋯ ∣ u n − 2 , n − 2 u n − 2 , n − 1 u n − 2 , n u n − 1 , n − 1 u n − 1 , n u n , n ] [ x 0 x 1 ∣ x n ] = [ y 0 y 1 ∣ y n ] x n = y n / u n n x n − 1 = ( y n − 1 − u n − 1 , n − 1 x n ) / u n − 1 , n − 1 ⋯ \begin{bmatrix} \smallsetminus & \cdots & \cdots& \lvert\\ & u_{n-2,n-2} & u_{n-2,n-1} & u_{n-2,n} \\ & & u_{n-1,n-1} & u_{n-1,n} \\ & & & u_{n,n} \end{bmatrix} \begin{bmatrix} x_{0}\\ x_{1}\\ \lvert\\ x_{n} \end{bmatrix} = \begin{bmatrix} y_{0}\\ y_{1}\\ \lvert\\ y_{n} \end{bmatrix}\\ \space \\ \space \begin{aligned} x_{n} &= y_{n} / u_{nn} \\ x_{n-1} &= (y_{n-1}-u_{n-1,n-1}x_{n}) / u_{n-1,n-1}\\ & \cdots \end{aligned} ⎣ ⎡∖⋯un−2,n−2⋯un−2,n−1un−1,n−1∣un−2,nun−1,nun,n⎦ ⎤⎣ ⎡x0x1∣xn⎦ ⎤=⎣ ⎡y0y1∣yn⎦ ⎤ xnxn−1=yn/unn=(yn−1−un−1,n−1xn)/un−1,n−1⋯

这样我们就得到了 x \bold x x的值。

2.1.1 直接法的优劣

对于比较稀疏的 A \bold A A， L \bold L L和 U \bold U U并不会太稀疏。 L \bold L L和 U \bold U U稀疏性取决于x元素的排列顺序（permutation）。

对于某些固定的 A \bold A A， L U \bf LU LU分解的计算和求解 x \bold x x的计算是分开的，因此分解后的 L \bold L L和 U \bold U U可以存起来复用。

算法在求解 x \bold x x时基本是无法并行的，因此不能用GPU加速。

2.2 迭代法

迭代法通常呈现如下的形式：

x [ k + 1 ] = x [ k ] + α M − 1 ( b − A x k ) \bold x^{[k+1]} = \bold x^{[k]} + \alpha \bold M^{-1}(\bold b - \bold A \bold x^{k}) x[k+1]=x[k]+αM−1(b−Axk)

α \alpha α称为松弛变量（Relaxation）。
M \bold M M称为迭代矩阵（Iterative Matrix）。
( b − A x k ) (\bold b - \bold A \bold x^{k}) (b−Axk)称为残差（Residual Error）。

算法的逻辑就是每次迭代都更新一次 x \bold x x， k k k就是迭代的次序。

理论上，如果求到了最终的解，那么残差项为 0 \bold 0 0。

2.2.1 收敛性证明

算法可以收敛的证明过程具体可以看课程21分46秒。

为了证明算法收敛，我们需要计算 b − A x k + 1 \bold b - \bold A \bold x^{k+1} b−Axk+1。

把上述的 x k + 1 \bold x^{k+1} xk+1的公式带入到 b − A x k + 1 \bold b - \bold A \bold x^{k+1} b−Axk+1中，得到：

b − A x k + 1 = b − A x k − α M − 1 ( b − A x k ) = ( I − α A M − 1 ) ( b − A x k ) \begin{aligned} \bold b - \bold A \bold x^{k+1} &= \bold b - \bold A \bold x^{k} - \alpha \bold M^{-1}(\bold b - \bold A \bold x^{k}) \\ &= (\bold I - \alpha \bold A\bold M^{-1})(\bold b - \bold A \bold x^{k}) \end{aligned} b−Axk+1=b−Axk−αM−1(b−Axk)=(I−αAM−1)(b−Axk)

再把 x k \bold x^{k} xk和 x k − 1 \bold x^{k-1} xk−1的关系公式带入，如此继续，最终得到

b − A x k + 1 = ( I − α A M − 1 ) k + 1 ( b − A x 0 ) \bold b - \bold A \bold x^{k+1} = (\bold I - \alpha \bold A\bold M^{-1})^{k+1}(\bold b - \bold A \bold x^{0}) b−Axk+1=(I−αAM−1)k+1(b−Ax0)

由于 b − A x 0 \bold b - \bold A \bold x^{0} b−Ax0是一个定值（因为算法一开始我们会给 x 0 \bold x^{0} x0赋一个初始值），所以只要证明 ( I − α A M − 1 ) k + 1 (\bold I - \alpha \bold A\bold M^{-1})^{k+1} (I−αAM−1)k+1无限趋近于0即可。

这里引入谱半径（Spectral Radius）的概念：

一个矩阵的谱半径，定义为某矩阵最大的特征值（eigen value）的绝对值，记作 ρ ( M ) \rho (\bold M) ρ(M)。

存在定理：

若矩阵 M \bold M M的 ρ ( M ) < 1 \rho (\bold M) < 1 ρ(M)<1，则 lim ⁡ k → + ∞ M k → 0 \lim\limits_{k \rarr +\infty} \bold M^{k} \rarr 0 k→+∞limMk→0。

综上所述，只要 ρ ( I − α A M − 1 ) < 1 \rho (\bold I - \alpha \bold A\bold M^{-1}) < 1 ρ(I−αAM−1)<1， b − A x k + 1 → 0 \bold b - \bold A \bold x^{k+1} \rarr 0 b−Axk+1→0，这样我们就可以得到正确的解。

所有的东西都归结于，我们如何构造迭代矩阵 M \bold M M上。

2.2.2 关于迭代矩阵

若M取作A的对角矩阵，该方法即为Jacobi Method，贾可比迭代。

若M去做A的下三角矩阵，该方法即为Gauss-Seidel Method。

还有切比雪夫（Chebyshev），共轭梯度（Conjugate Gradient）等方法。

2.2.3 迭代法的优劣

优点

简便
对于不是十分精确的求解过程来说，快
可并行

缺点

有收敛条件
精确解最好不使用

【GAMES103学习笔记】线性系统（Linear System）及其解法（Linear Sovler）相关推荐

sv_labs学习笔记——sv_lab5_下(System Verilog)
本文延续前一篇sv_labs学习笔记--sv_lab5_上(System Verilog),进一步学习完善lab5的内容 sv_labs学习笔记--sv_lab5_下(System Verilog) ...
sv_labs学习笔记——sv_lab5_上(System Verilog)
本节将介绍lab5的第一部分,主要总结一般设计学习与思考的方式与需要着重学习的点,同时以lab5作为参考,分析数据流流向,验证组件的通信与抽象化,实现的整体思路. sv_labs学习笔记--sv_la ...
【GAMES103学习笔记】刚体（Rigid Body）
[GAMES103学习笔记]刚体(Rigid Body) 0 什么是刚体及刚体模拟 1 平移 1.1 平移运动 1.1.1 速度 1.1.2 位置 1.2 时间积分 1.3 半隐式欧拉 1.4 力的分 ...
C#学习笔记：利用System,EventArgs实现委托，响应键盘按键事件
参考书目:C#6.0学习笔记--从第一行C#代码到第一个项目设计(作者周家安)P96 学习目的:掌握System,EventArgs实现委托的方法,响应键盘按键事件.捕捉用户的键盘输入,然后触发Key ...
FPGA学习笔记_Quartus II_In system sources and probes editor(ISSP)调试工具的使用
FPGA学习笔记 Quartus II prime Standard Edition-In system sources and probes editor(ISSP)调试工具的使用 Quartus ...
Java基础学习笔记之：System类;Math类;Arrays类BigInteger，BigDecimal
System类在API中System类介绍的比较简单,我们给出定义,System中代表程序所在系统,提供了对应的一些系统属性信息,和系统操作. System类不能手动创建对象,因为构造方法被priv ...
部份API学习笔记（Math,System,Object,Date,SimpleDateFormat）
一.Math类:Math包含执行基本数字运算的方法,如基本指数,对数,平方根和三角函数.是静态类,用static修饰的,没有构造方法,不有实例化对象,直接用类名调用方法例属于:java.lang 使用 ...
算法学习笔记2：凸包及其解法
凸包及其解法凸包概念解法暴力解法判断方向分治点到直线距离 Jarvis步进法 Graham扫描法 Andrew算法凸包概念在一个实数向量空间V中,对于给定集合X,所有包含X的凸集的交集 ...
APUE 学习笔记 - Chapter 6. System Data File and Infomation
1.密码文件每个系统都会有一个文件统一记录用户名与密码,通常是/etc/passwd. 关于这个文件有: root 的 uin 通常为 0 . 文 ...

【GAMES103学习笔记】线性系统（Linear System）及其解法（Linear Sovler）

【GAMES103学习笔记】线性系统（Linear System）及其解法（Linear Sovler）

0. 为什么要研究线性系统

1. 什么是线性系统

2. 解线性系统

2.1 直接法

2.1.1 直接法的优劣

2.2 迭代法

2.2.1 收敛性证明

2.2.2 关于迭代矩阵

2.2.3 迭代法的优劣

【GAMES103学习笔记】线性系统（Linear System）及其解法（Linear Sovler）相关推荐

最新文章

热门文章