mirror of
https://github.com/cheetahlou/CategoryResourceRepost.git
synced 2025-11-17 22:53:42 +08:00
mod
This commit is contained in:
louzefeng
200
极客时间专栏/重学线性代数/应用篇/11 | 如何运用线性代数方法解决图论问题？.md
Normal file
200
极客时间专栏/重学线性代数/应用篇/11 | 如何运用线性代数方法解决图论问题？.md
Normal file
@@ -0,0 +1,200 @@
|
||||
<audio id="audio" title="11 | 如何运用线性代数方法解决图论问题？" controls="" preload="none"><source id="mp3" src="https://static001.geekbang.org/resource/audio/a5/9a/a56b3e2600f6e91eff82a1c269df3a9a.mp3"></audio>
|
||||
|
||||
你好，我是朱维刚。欢迎你继续跟我学习线性代数，今天我要讲的内容是“如何运用线性代数方法解决图论问题”。
|
||||
|
||||
“图”这个字在计算机科学领域很常见，特别是在数据结构中。一说到图，是必定要联系到图论（Graph Theory）的，因为它是以图为研究对象的数学的一个分支。图论中的图，是由若干给定的**点**及连接两点的**线**所构成的图形，这种图形通常用来描述某些事物之间的某种特定关系，用点代表事物，用连接两点的线表示相应两个事物间具有这种关系。
|
||||
|
||||
说到这，你也许会问，这个和线性代数、矩阵有什么关系？
|
||||
|
||||
## 图的数学定义
|
||||
|
||||
既然是数学课，我们还是要先讲一下图的数学定义：一个图$G$是指一个有序三元组$(V, E, \phi)$，$V$是非空的顶点集；$E$是不与$V$相交的边集；$\phi$是关联函数，它使$G$的每条边对应于$G$的无序顶点对。如果$e$是一条边，$u$和$v$是顶点，使得$\phi(e)=u v$，则$e$连接$u$和$v$，也就是顶点$u$和$v$是$e$的端点。
|
||||
|
||||
好了，现在是时候通过两个应用场景来看下，如何把矩阵和图论关联起来，并运用在解决实际问题中了。
|
||||
|
||||
## 邻接矩阵应用
|
||||
|
||||
首先，是邻接矩阵（Adjacency Matrix），邻接矩阵是表示顶点之间相邻关系的矩阵。假设$G$是一个图，$V(G)$是$G$的顶点集，$E(G)$是$G$的边集，设$G$中有$n$个顶点，$v_{1}, v_{2}, \ldots, v_{n}$，$A=(a_{ij})_{n×n}$是$G$的邻接矩阵。
|
||||
|
||||
$$a_{i}=\left\{\begin{array}{l}<br>
|
||||
1, v_{i} v_{j} \in E(G) \\\<br>
|
||||
0, v_{i} v_{j} \notin E(G)^{\prime}, i, j=1,2, \ldots, n<br>
|
||||
\end{array}\right.$$
|
||||
|
||||
已知情况是这样的，那我们现在来看一下邻接矩阵在现实问题中的应用。这个例子来源于1994年全国大学生数学建模竞赛试题B题。
|
||||
|
||||
某厂生产一种弹子锁具，每个锁具的钥匙有5个槽，每个槽的高度从$｛1,2,3,4,5,6｝$中任取一数。由于工艺及其他原因，制造锁具时对5个槽的高度还有两个限制。
|
||||
|
||||
1. 至少有3个不同的数；
|
||||
1. 相邻两槽的高度之差不能为5。
|
||||
|
||||
满足以上条件制造出来的所有互不相同的锁具称为一批。销售部门在一批锁具中随意地取每60个装成一箱出售。问：每一批锁具有多少个，装多少箱？
|
||||
|
||||
我们先来看下弹子锁具的样子，否则自己想象会要一些时间。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/1b/04/1b45eea656a8cc18dde70c16f5a17204.gif" alt="">
|
||||
|
||||
这是一个7个槽的弹子锁具，只是比我们例子的5个槽多了两个槽。有了弹子锁具形象化输出后，我们开始解题。锁具装箱是一个排列组合的数学问题，但如果我们用图论的邻接矩阵方法来解这个问题，就能够起到化繁为简的作用。
|
||||
|
||||
首先，我们构造一个6节点的图：把1、2、3、4、5、6这6个数作为6个节点，如果两个数字可以相邻，那这两个节点之间就加一条边，每个节点有自己到自己的一条边。于是，我们得到了锁具各槽之间的关系图：
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/ba/0f/ba703a93fec043yyf741f221e08ca50f.png" alt="">
|
||||
|
||||
接着，构建邻接矩阵$A$，根据前面说的，如果两点之间有一条边，那在矩阵中，相应位置的值就是1，比如：节点1和2之间有一条边，那矩阵第一行第二列和第二行第一列的值就是1，节点1和6之间没有边，那矩阵第一行第六列和第六行第一列的值就是0，因为每个节点有自己到自己的一条边，所以第一行第一列的值就是1，其它5个节点也是一样的。
|
||||
|
||||
$$A=\left[\begin{array}{llllll}<br>
|
||||
1 & 1 & 1 & 1 & 1 & 0 \\\<br>
|
||||
1 & 1 & 1 & 1 & 1 & 1 \\\<br>
|
||||
1 & 1 & 1 & 1 & 1 & 1 \\\<br>
|
||||
1 & 1 & 1 & 1 & 1 & 1 \\\<br>
|
||||
1 & 1 & 1 & 1 & 1 & 1 \\\<br>
|
||||
0 & 1 & 1 & 1 & 1 & 1<br>
|
||||
\end{array}\right]$$
|
||||
|
||||
因为我们从没有1、6相邻的关系图得到了邻接矩阵$A$，所以$A$中所有元素之和表示两个槽高**无1、6相邻**的锁具个数。而每个无1、6相邻的5位数与关系图中长度是1的一条链一一对应。于是，$A^{k}$中各元素之和就是长度为$k$的链接个数。比如，$A^{2}$中第$i$行第$j$列的元素就是$i$开始经过两条边到达$j$的链接个数。我们这里因为是5个元素，也就是要经过4条边，所以需要计算$A^{4}$。
|
||||
|
||||
$$A^{4}=\left[\begin{array}{cccccc}<br>
|
||||
141 & 165 & 165 & 165 & 165 & 140 \\\<br>
|
||||
165 & 194 & 194 & 194 & 194 & 165 \\\<br>
|
||||
165 & 194 & 194 & 194 & 194 & 165 \\\<br>
|
||||
165 & 194 & 194 & 194 & 194 & 165 \\\<br>
|
||||
165 & 194 & 194 & 194 & 194 & 165 \\\<br>
|
||||
140 & 165 & 165 & 165 & 165 & 141<br>
|
||||
\end{array}\right]$$
|
||||
|
||||
把$A^{4}$中的元素求和，就能得到相邻高差为5的锁具数是6306。
|
||||
|
||||
最后，因为题目的限制提到了槽的高度至少有3个不同的数，我们还要把6306这个数字减去仅有一个、两个槽高的锁具：$6306-\left(6+\left(C_{6}^{2}-1\right)\left(2^{5}-2\right)\right)=5880$。
|
||||
|
||||
所以，我们得到一批锁具的个数是5880，总共装5880/60=98箱。
|
||||
|
||||
这样，我们通过邻接矩阵的图论知识，解决了一批锁具的数量问题，比其它方法看起来更简单。
|
||||
|
||||
>
|
||||
特别提示：文中用到的$A^{k}$在图论中的实际意义，是来自刘亚国的一篇文献《图论中邻接矩阵的应用》。
|
||||
|
||||
|
||||
## 关联矩阵应用
|
||||
|
||||
接下来，我们在邻接矩阵上再升级一下，把边变成有向边。这样就形成了另一类矩阵：关联矩阵。关联矩阵经常被用在图论中，现在我们就来看一下关联矩阵和图之间的关系。如下图所示，我们定义了一个拥有4个节点和6条边的图。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/d4/f6/d447367b8b2980d05f44aeac421664f6.png" alt="">
|
||||
|
||||
接下来，定义一个6×4的矩阵来描述这个图，其中列表示点(1)到点(4)，行表示边1到边6：
|
||||
|
||||
$$A=\left[\begin{array}{cccc}<br>
|
||||
-1 & 1 & 0 & 0 \\\<br>
|
||||
-1 & 0 & 1 & 0 \\\<br>
|
||||
0 & -1 & 1 & 0 \\\<br>
|
||||
-1 & 0 & 0 & 1 \\\<br>
|
||||
0 & -1 & 0 & 1 \\\<br>
|
||||
0 & 0 & -1 & 1<br>
|
||||
\end{array}\right]$$
|
||||
|
||||
矩阵$A$只包含了三类元素：-1、1、0。-1表示点的箭头的出方向，1表示点的箭头的入方向，0则表示点和点之间没有关联。举例来说，矩阵$A$的第一行元素是-1、1、0、0，那对于边1和点(1)、点(2)说，我们从图中可以看到边1是从点(1)到点(2)，$A$中-1对于点(1)来说是箭头的出方向，1对于点(2)来说是箭头的入方向，而边1和点(3)、点(4)没有任何关系，所以第一行第三列和第四列都是0。
|
||||
|
||||
这里，我们把关联矩阵用到现实场景中，比如：让它为电子电路服务，用它来分析整个电路的情况，也就是电路的拓扑结构，这里的电路指的是基尔霍夫定律，是分析和计算较为复杂电路的基础。假设$x_{1}, x_{2}, x_{3}, x_{4}$是这几个点的电压值，我们来看一下$Ax$的结果：
|
||||
|
||||
$$A x=\left[\begin{array}{cccc}<br>
|
||||
-1 & 1 & 0 & 0 \\\<br>
|
||||
-1 & 0 & 1 & 0 \\\<br>
|
||||
0 & -1 & 1 & 0 \\\<br>
|
||||
-1 & 0 & 0 & 1 \\\<br>
|
||||
0 & -1 & 0 & 1 \\\<br>
|
||||
0 & 0 & -1 & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{l}<br>
|
||||
x_{1} \\\<br>
|
||||
x_{2} \\\<br>
|
||||
x_{3} \\\<br>
|
||||
x_{4}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
x_{2}-x_{1} \\\<br>
|
||||
x_{3}-x_{1} \\\<br>
|
||||
x_{3}-x_{2} \\\<br>
|
||||
x_{4}-x_{1} \\\<br>
|
||||
x_{4}-x_{2} \\\<br>
|
||||
x_{4}-x_{3}<br>
|
||||
\end{array}\right]$$
|
||||
|
||||
由此可见，结果是差值，也就是沿着边1到6的电势差，有了电势差，就说明有电流，但如果Ax=0会怎样呢？也就是方程满足这样的等式：
|
||||
|
||||
$$A x=\left[\begin{array}{cccc}<br>
|
||||
-1 & 1 & 0 & 0 \\\<br>
|
||||
-1 & 0 & 1 & 0 \\\<br>
|
||||
0 & -1 & 1 & 0 \\\<br>
|
||||
-1 & 0 & 0 & 1 \\\<br>
|
||||
0 & -1 & 0 & 1 \\\<br>
|
||||
0 & 0 & -1 & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{l}<br>
|
||||
x_{1} \\\<br>
|
||||
x_{2} \\\<br>
|
||||
x_{3} \\\<br>
|
||||
x_{4}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
x_{2}-x_{1} \\\<br>
|
||||
x_{3}-x_{1} \\\<br>
|
||||
x_{3}-x_{2} \\\<br>
|
||||
x_{4}-x_{1} \\\<br>
|
||||
x_{4}-x_{2} \\\<br>
|
||||
x_{4}-x_{3}<br>
|
||||
\end{array}\right]=\left[\begin{array}{l}<br>
|
||||
0 \\\<br>
|
||||
0 \\\<br>
|
||||
0 \\\<br>
|
||||
0 \\\<br>
|
||||
0 \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]$$
|
||||
|
||||
很明显，这些差值，也就是电势差都等于0，意味着没有电流。同理，如果把电压值换成温度呢？那应用场景就切换成热传导了。
|
||||
|
||||
刚才我们看到了$Ax=0$的情况，你还记得[第八篇](https://time.geekbang.org/column/article/272815)中说的零空间吗？它关注的就是$Ax=0$，也就是向量空间$V$中所有经过$\phi$映射为零的向量集合，用符号表示就是：$ker(\phi)$，它的维数叫做零化度（nullity），表示成：$dim(ker(\phi))$。
|
||||
|
||||
而在电路例子中，它表示的是所有六个电势差都是0，也就意味着：所有四个电压值是相等的，在零空间中的每个$x$都是一个常向量：$x=(c,c,c,c)$。所以，$A$的零空间是一条线。无论我们怎么同时升高或降低电压量$c$，都不会改变电势差0。
|
||||
|
||||
刚才我们说的是电压，现在我们来具体看看关联矩阵在基尔霍夫电流定律上的运用。我们来定义一个拥有4个节点和5条边的图：
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/ce/d6/ce37e3d70213cfc23f6c33eb00637ed6.png" alt="">
|
||||
|
||||
基尔霍夫电流定律定义：$A^{T} y=0$，其中y是向量$y_{1}, y_{2}, y_{3}, y_{4}, y_{5}$，我们把这个图以关联矩阵的形式写出来就是：
|
||||
|
||||
$$\left[\begin{array}{ccccc}<br>
|
||||
-1 & 0 & -1 & -1 & 0 \\\<br>
|
||||
1 & -1 & 0 & 0 & 0 \\\<br>
|
||||
0 & 1 & 1 & 0 & -1 \\\<br>
|
||||
0 & 0 & 0 & 1 & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{l}<br>
|
||||
y_{1} \\\<br>
|
||||
y_{2} \\\<br>
|
||||
y_{3} \\\<br>
|
||||
y_{4} \\\<br>
|
||||
y_{5}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
0 \\\<br>
|
||||
0 \\\<br>
|
||||
0 \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]$$
|
||||
|
||||
这里-1，0，1的含义上面有所描述，第一行和$y$向量相乘后得到：$-y_{1}-y_{3}-y_{4}=0$，说明从节点1出来的总电流等于0，满足守恒定律；第二行和$y$向量相乘后得到：$y_{1}-y_{2}=0$，说明流入节点2的电流和从节点2流出的电流相等；同样，后面两行分别和$y$向量相乘后得到：$y_{2}+y_{3}-y_{5}=0$和$y_{4}+y_{5}=0$，和图表示的都一致，也都符合守恒定律。
|
||||
|
||||
好了，到这里简单电路的数学知识，也就是关联矩阵讲完了，如果碰到更复杂的电路，比如：在节点之间，也就是边上有电流源，那么，等式就要从$A^{T} y=0$变成$A^{T} y=f$。
|
||||
|
||||
## 本节小结
|
||||
|
||||
本节是第一篇应用篇，所以我从更贴近生活的例子来讲解线性代数的应用，通过弹子锁具，让你能够了解，邻接矩阵与图论之间是怎么关联的；通过基尔霍夫定律，让你能够了解关联矩阵与图论之间是怎么关联的。
|
||||
|
||||
所以，邻接矩阵、关联矩阵的最大意义就是用数学的方式描述图，进而来描述某些事物之间的某种特定关系，是不是发现问题后通过数学工具来解决问题很美妙呢？
|
||||
|
||||
## 线性代数练习场
|
||||
|
||||
这次的练习题稍微有些难度，是一道传统练习题。
|
||||
|
||||
三名商人各带一个随从乘船渡河，现有一只小船只能容纳两个人，由他们自己划行，若在河的任一岸的随从人数多于商人，他们就可能抢劫财物。但如何乘船渡河由商人决定，试给出一个商人安全渡河的方案，使得渡河的次数最少。
|
||||
|
||||
>
|
||||
<p>注意：这里的问题包含两层含义——安全渡河和渡河次数最少。<br>
|
||||
提示：使用本节的第一个例子的邻接矩阵和$A^{k}$来解这道题。</p>
|
||||
|
||||
|
||||
欢迎在留言区晒出你的运算过程和结果，我会及时回复。同时，也欢迎你把这篇文章分享给你的朋友，一起讨论、学习。
|
||||
247
极客时间专栏/重学线性代数/应用篇/12 | 如何通过矩阵转换让3D图形显示到二维屏幕上？.md
Normal file
247
极客时间专栏/重学线性代数/应用篇/12 | 如何通过矩阵转换让3D图形显示到二维屏幕上？.md
Normal file
@@ -0,0 +1,247 @@
|
||||
<audio id="audio" title="12 | 如何通过矩阵转换让3D图形显示到二维屏幕上？" controls="" preload="none"><source id="mp3" src="https://static001.geekbang.org/resource/audio/f4/1f/f41d99f188c0c728419b92967449da1f.mp3"></audio>
|
||||
|
||||
你好，我是朱维刚。欢迎你继续跟我学习线性代数，今天我要讲的内容是“如何通过矩阵转换让3D图形显示到二维屏幕上”。
|
||||
|
||||
在第八篇的[线性映射](https://time.geekbang.org/column/article/272815)中，我从二维直角坐标系的角度，讲解了线性映射和变换矩阵。其中，我特别讲到了，二维平面图形图像处理中的线性变换，比如物体的拉伸和旋转。在第九篇的[仿射空间](https://time.geekbang.org/column/article/274222)中，更是提到了3D的平移矩阵、缩放矩阵和旋转矩阵。
|
||||
|
||||
而这一篇则有些不一样，我会从更实践的角度，让你了解到二维平面和三维空间的变换，以及3D图形是如何显示到二维屏幕上的。矩阵在这里扮演的角色可以说是功不可没，接下来我们一起来看下矩阵到底是怎么做到的。
|
||||
|
||||
## 三维空间变换
|
||||
|
||||
我们都知道，计算机图形图像处理的是图片，且计算机屏幕是二维的。那你有没有想过，我们在屏幕上看到的静态和动态三维世界到底是怎么回事呢？这个就要涉及到三维到二维的投影技术了，这类技术都离不开矩阵，而且是超大规模矩阵运算。
|
||||
|
||||
三维空间的变换依赖于4×4矩阵，可能你会想，为什么不是3×3呢？这是因为四个关键运算中有一个无法用3×3矩阵来完成，其他三个运算为了统一也就都采用4×4矩阵了，这四个关键运算是：
|
||||
|
||||
- 平移；
|
||||
- 缩放；
|
||||
- 旋转；
|
||||
- 投影。
|
||||
|
||||
平移就是那个无法用3×3矩阵来完成的特殊运算，也是看起来最简单的运算，只是每个点都加上向量$v_{0}$，也就是点$(x_{0},y_{0},z_{0})$。
|
||||
|
||||
但是，你别被这个假象欺骗了，平移这个运算是非线性的。这一点只需要看平移前各点与原点的连线，以及平移后各点与原点之间的连线就知道了。或者，你也可以从公式的角度理解，就是$f(a+b)$不等于$f(a)+f(b)$。而为了表示平移，以及现实世界的描述，就需要使用第九篇中说的[仿射空间](https://time.geekbang.org/column/article/274222)。所以，3×3矩阵是无法平移原点的。
|
||||
|
||||
但是，如果我们把原点坐标变成$(0,0,0,1)$，那就能解决平移的问题了。点$(x,y,z)$的齐次坐标就是$(x,y,z,1)$，这就变成了4×4矩阵。接下来，我分别介绍这四个关键运算，它们是3D图形显示在屏幕上的第一步，也就是坐标系变换要做的事情，比如：将一个点从局部坐标系变换到世界坐标系是通过平移、缩放及旋转矩阵进行的。
|
||||
|
||||
## 平移
|
||||
|
||||
我们沿着向量$v_{0}$平移整个三维空间，把原点平移到了$(x_{0},y_{0},z_{0})$，这也就意味着三维空间的每个点都加上了点$(x_{0},y_{0},z_{0})$。使用齐次坐标，把整个空间平移了$v_{0}$的4×4矩阵$T$如下所示。
|
||||
|
||||
$$<br>
|
||||
T=\left[\begin{array}{llll}<br>
|
||||
1 & 0 & 0 & 0 \\\<br>
|
||||
0 & 1 & 0 & 0 \\\<br>
|
||||
0 & 0 & 1 & 0 \\\<br>
|
||||
x_{0} & y_{0} & z_{0} & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
这里很重要的一点是，计算机图形图像是基于行向量计算的。也就是说，计算方法是行乘矩阵，而不是矩阵乘列，比如：$\left[\begin{array}{llllllll}0 & 0 & 0 & 1\end{array}\right] T=\left[\begin{array}{llll}x_{0} & y_{0} & z_{0} & 1\end{array}\right]$。
|
||||
|
||||
平移的整个过程是这样的：假设要把原来的某个点$(x,y,z)$平移$v_{0}$，我们需要切换到齐次坐标$(x,y,z,1)$，然后，$(x,y,z,1)$再乘$T$，就能得到每个原来的向量$v$平移到$v+v_{0}$的最终结果：$\left[\begin{array}{llll}x & y & z & 1\end{array}\right] T=\left[\begin{array}{lllll}x+x_{0} & y+y_{0} & z+z_{0} & 1\end{array}\right]$。
|
||||
|
||||
这里你需要注意：一个行向量乘T的结果还是一个行向量。
|
||||
|
||||
## 缩放
|
||||
|
||||
在前端开发中，我们经常会调整图片宽度和高度来适配页面，比如：把图片整体放大90%，那么在线性代数中就是0.9乘单位矩阵。在二维平面中，我们通常用2×2矩阵来表达缩放，在三维立体中则是3×3矩阵。而在计算机图形图像的齐次坐标中，就不一样了，需要大一个维度，也就是说，3×3矩阵变成了4×4矩阵。
|
||||
|
||||
比如，二维平面中图片放大90%就是：
|
||||
|
||||
$$<br>
|
||||
S=\left[\begin{array}{ccc}<br>
|
||||
0.9 & 0 & 0 \\\<br>
|
||||
0 & 0.9 & 0 \\\<br>
|
||||
0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
三维立体中图片放大90%就是：
|
||||
|
||||
$$<br>
|
||||
S=\left[\begin{array}{cccc}<br>
|
||||
0.9 & 0 & 0 & 0 \\\<br>
|
||||
0 & 0.9 & 0 & 0 \\\<br>
|
||||
0 & 0 & 0.9 & 0 \\\<br>
|
||||
0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
缩放还可以在不同的方向上进行，比如：一个二维平面图片从整页适配调整到半页适配，$y$方向就要乘$\frac{1}{2}$，创建一个$\frac{1}{4}$的页边留白，$x$方向就要乘$\frac{3}{4}$，这样得到的缩放矩阵就是：
|
||||
|
||||
$$<br>
|
||||
S=\left[\begin{array}{lll}<br>
|
||||
\frac{3}{4} & 0 & 0 \\\<br>
|
||||
0 & \frac{1}{2} & 0 \\\<br>
|
||||
0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
平移和缩放组合情况会怎样呢？如果我们要先平移再缩放，那应该这样乘：$vTS$，如果我们要先缩放再平移，那应该这样乘：$vST$。注意：它们乘的顺序是不同的，哪个运算先做就先乘，因为矩阵的左乘和右乘的结果是不同的。
|
||||
|
||||
在第九篇的[仿射空间](https://time.geekbang.org/column/article/274222)中提到了平移和缩放矩阵，你也可以回过头再去看看。
|
||||
|
||||
## 旋转
|
||||
|
||||
二维和三维空间的旋转由正交矩阵$Q$来完成，它的行列式是+1。同样我们使用齐次坐标，一个平面旋转的正交矩阵$Q$就从2×2就变成了3×3矩阵$R$。
|
||||
|
||||
$$<br>
|
||||
Q=\left[\begin{array}{cc}<br>
|
||||
\cos \theta & -\sin \theta \\\<br>
|
||||
\sin \theta & \cos \theta<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
R=\left[\begin{array}{ccc}<br>
|
||||
\cos \theta & -\sin \theta & 0 \\\<br>
|
||||
\sin \theta & \cos \theta & 0 \\\<br>
|
||||
0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
这个矩阵是围绕原点旋转了平面，那如果矩阵旋转时围绕的不是原点，而是其他点呢？这个就稍微复杂一些，不是直接旋转，而是先平移再旋转，比如我们要围绕点$(4,5)$，让平面旋转$\theta$角度的话：
|
||||
|
||||
1. 首先，要把$(4,5)$平移到$(0,0)$；
|
||||
1. 接着，旋转$\theta$角度；
|
||||
1. 最后，再把$(0,0)$平移回$(4,5)$。
|
||||
|
||||
整个过程通过数学公式来表达就是：
|
||||
|
||||
$$<br>
|
||||
v T_{00} R T_{45}=\left[\begin{array}{lll}<br>
|
||||
x & y & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{ccc}<br>
|
||||
1 & 0 & 0 \\\<br>
|
||||
0 & 1 & 0 \\\<br>
|
||||
-4 & -5 & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{ccc}<br>
|
||||
\cos \theta & -\sin \theta & 0 \\\<br>
|
||||
\sin \theta & \cos \theta & 0 \\\<br>
|
||||
0 & 0 & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{ccc}<br>
|
||||
1 & 0 & 0 \\\<br>
|
||||
0 & 1 & 0 \\\<br>
|
||||
4 & 5 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
说完二维我们再来说三维。不过在三维空间中，旋转就有些不一样了，因为它是围绕一个轴“翻转”的。更“数学”的说法就是，围绕$λ=1$的特征向量的一条线翻转。
|
||||
|
||||
现在，我们来看看分别围绕$x$、$y$和$z$轴方向旋转的矩阵$R$有什么不同？
|
||||
|
||||
1.围绕$x$轴方向旋转：
|
||||
|
||||
$$<br>
|
||||
R_{x}=\left[\begin{array}{cccc}<br>
|
||||
1 & 0 & 0 & 0 \\\<br>
|
||||
0 & \cos \theta & -\sin \theta & 0 \\\<br>
|
||||
0 & \sin \theta & \cos \theta & 0 \\\<br>
|
||||
0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
2.围绕$y$轴方向旋转：
|
||||
|
||||
$$<br>
|
||||
R_{y}=\left[\begin{array}{cccc}<br>
|
||||
\cos \theta & 0 & \sin \theta & 0 \\\<br>
|
||||
0 & 1 & 0 & 0 \\\<br>
|
||||
-\sin \theta & 0 & \cos \theta & 0 \\\<br>
|
||||
0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
3.围绕$z$轴方向旋转：
|
||||
|
||||
$$<br>
|
||||
R_{z}=\left[\begin{array}{cccc}<br>
|
||||
\cos \theta & -\sin \theta & 0 & 0 \\\<br>
|
||||
\sin \theta & \cos \theta & 0 & 0 \\\<br>
|
||||
0 & 0 & 1 & 0 \\\<br>
|
||||
0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
你看出来哪里不同了吗？其实主要就是1的位置不同，以及$y$轴方向旋转的$sin$互换了。
|
||||
|
||||
## 投影
|
||||
|
||||
现在，我们想把3D图形显示到二维屏幕上，该怎么做呢？
|
||||
|
||||
从数学角度理解就是把三维向量投影到平面上。在线性代数中，我们看到的大部分的平面都是通过原点的，但在现实生活中则不是。一个通过原点的平面是一个向量空间，而其他的平面则是仿射空间，具体仿射空间的定义你可以回顾一下[第九篇](https://time.geekbang.org/column/article/274222)的内容。
|
||||
|
||||
我们先来看看平面通过原点的情况。假设一个通过原点的平面，它的单位法向量是$n$，那么平面中的向量$v$，满足这个等式：$n^{T}v=0$。
|
||||
|
||||
而投影到平面的投影矩阵是：$I-nn^{T}$。
|
||||
|
||||
如果把原来的向量和这个投影矩阵相乘，就能投影这个向量。我们可以用这个投影矩阵来验证一下：单位法向量$n$投影后成为了0向量，而平面向量$v$投影后还是其自身。
|
||||
|
||||
$$<br>
|
||||
(I-n n^{T}) n=n-n(n^{T} n)=0<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
(I-n n^{T}) v=v-n(n^{T} v)=v<br>
|
||||
$$
|
||||
|
||||
接下来，我们在齐次坐标中来看一下4×4的投影矩阵：
|
||||
|
||||
$$<br>
|
||||
P=\left[\begin{array}{lll}<br>
|
||||
& & & 0 \\\<br>
|
||||
& I-n n^{T} & & 0 \\\<br>
|
||||
& & & 0 \\\<br>
|
||||
0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
假设现在有一个不过原点的平面，$v_{0}$是这个平面上的一个点，现在要把$v_{0}$投影到这个平面，则需要经历三个步骤，和刚才介绍的围绕点$(4,5)$，让平面旋转$\theta$角度经历的三个步骤类似：
|
||||
|
||||
1. 把$v_{0}$平移到原点；
|
||||
1. 沿着$n$方向投影；
|
||||
1. 再平移回$v_{0}$。
|
||||
|
||||
整个过程通过数学公式来表达就是：
|
||||
|
||||
$$<br>
|
||||
T_{-v_{0}} P T_{+v_{0}}=\left[\begin{array}{cc}<br>
|
||||
I & 0 \\\<br>
|
||||
-v_{0} & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{cc}<br>
|
||||
I-n n^{T} & 0 \\\<br>
|
||||
0 & 1<br>
|
||||
\end{array}\right]\left[\begin{array}{ll}<br>
|
||||
I & 0 \\\<br>
|
||||
v_{0} & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
## 计算机3D图形介绍
|
||||
|
||||
有了数学知识的铺垫，我们再来看计算机3D图形显示到二维屏幕上的过程。在3D环境中，三维物体从取景到屏幕显示，需要经历一系列的坐标变换（又称为空间变换），才能生成二维图像显示在输出设备上。
|
||||
|
||||
将一个3D物体显示出来需要经历三个步骤，其中，第一步，也是最重要的一步就是坐标系变换，将局部坐标系表示的点变换到世界坐标系中，然后再变换到视图坐标系（或叫摄像机坐标系），接着继续变换到裁剪坐标系（投影坐标系）。
|
||||
|
||||
- 将一个点从局部坐标系变换到世界坐标系是通过平移、缩放及旋转矩阵进行的。
|
||||
- 如果将世界坐标系中的一个点变换到视图坐标系（摄像机坐标系），则可以使用视图矩阵进行操作。视图矩阵我们这里没有详细说明，它有个相对复杂的推导过程的，感兴趣的同学可以参考我后面推荐的两本书。
|
||||
- 如果将视图坐标系（摄像机坐标系）中的一个点变换到裁剪坐标系（投影坐标系），则可以使用投影矩阵进行操作。
|
||||
|
||||
最后，我推荐两本非常好的书作为你继续研究计算机3D图形的参考。
|
||||
|
||||
>
|
||||
<p>《TypeScript图形渲染实战：基于WebGL的3D架构与实现》，作者：步磊峰，这本书描述了3D图形处理的基本数学知识的同时，更注重WebGL框架下的图形渲染实战。<br>
|
||||
《Computer Graphics: Principles and Practice (3rd Edition)》，作者：Hughes, Van Dam, McGuire, Skylar, Foley, Feiner, Akeley，这本书虽然也有实践，但更偏重计算机图形理论一些。</p>
|
||||
|
||||
|
||||
## 本节小结
|
||||
|
||||
今天的整篇内容都是围绕三维空间的变换展开的，你需要掌握三维空间中的四个关键运算：平移、缩放、旋转和投影的基本概念，以及对应的平移、缩放、旋转和投影矩阵，这些都是继续深入学习计算机3D图形处理的数学基础。
|
||||
|
||||
因为在3D环境中，三维物体从取景到屏幕显示，需要经历一系列的坐标变换，才能生成二维图像显示在输出设备上。了解了这些之后，你就能掌握计算机3D图形处理的本质，也许还能在将来的实践中优化图形渲染效率。
|
||||
|
||||
## 线性代数练习场
|
||||
|
||||
今天我要给你一道开放题：如果把正方形投影到一个平面上，你会得到一个什么形状的图形？
|
||||
|
||||
欢迎在留言区晒出你的结果和思考，我会及时回复。同时，也欢迎你把这篇文章分享给你的朋友，一起讨论、学习。
|
||||
353
极客时间专栏/重学线性代数/应用篇/13 | 如何通过有限向量空间加持的希尔密码，提高密码被破译的难度？.md
Normal file
353
极客时间专栏/重学线性代数/应用篇/13 | 如何通过有限向量空间加持的希尔密码，提高密码被破译的难度？.md
Normal file
@@ -0,0 +1,353 @@
|
||||
<audio id="audio" title="13 | 如何通过有限向量空间加持的希尔密码，提高密码被破译的难度？" controls="" preload="none"><source id="mp3" src="https://static001.geekbang.org/resource/audio/6d/bb/6da29f4efc61fec95f9bdcc6915413bb.mp3"></audio>
|
||||
|
||||
你好，我是朱维刚。欢迎你继续跟我学习线性代数。
|
||||
|
||||
今天我要讲的内容是“如何通过有限向量空间加持的希尔密码，提高密码被破译的难度”。
|
||||
|
||||
这篇的内容会非常有趣，是和密码加密、解密有关的。不知道你有没有看过电影《模仿游戏》，故事描述的是阿兰·图灵在二战期间破译德军的恩尼格玛密码机（Enigma），很精彩，我看了很多遍。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/ac/5b/ac9afcc1194bc25f2ddcc5fb109bb85b.jpg" alt="">
|
||||
|
||||
不过电影毕竟是电影，有许多内容是不现实的，好在表达出来的破译恩尼格玛密码的核心观点是正确的。要破译一份被恩尼格玛机加密的密文，需要这三类信息：
|
||||
|
||||
1. 恩格玛机的工作原理及内部构造，包括每个转子的线路连接；
|
||||
1. 德军对恩格玛机的操作守则；
|
||||
1. 德军所使用的每日初始设置。恩格玛机的每日初始设置包含了三个信息：即转子的排列顺序、每个转子的初始位置，以及插线板的设置。这些信息被印刷在密码本上分发至德军全军，每24小时更换一次设置，每月更换一次密码本。
|
||||
|
||||
这些在电影里确实都交代了，我也不过多剧透了。其实，恩尼格玛密码机的本质就是**替换密码**。而今天我要讲的也是一种替换密码——希尔密码。因为我们专栏讲的是线性代数，所以，这篇应用我们会以矩阵论原理为基础，来进行讲解。
|
||||
|
||||
## 为什么需要希尔密码？
|
||||
|
||||
要讲密码，我们得先知道人们为什么需要它。
|
||||
|
||||
最古老、最原始的加密算法，会把明文的字母按照某种配对关系替换成其他的字母，从而得到一段别人看不懂的密文，许多谍战剧用到过这类方法。看起来，这个方法好像很难人为进行破解，但从语言和统计学角度看，它其实是漏洞百出的。
|
||||
|
||||
举个例子，在一篇普通英语文章中，各字母出现的概率有很大的不同。如果我们对足够多的文本进行分析，就可以统计出每一个字母在英文文本中出现的平均概率。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/6d/af/6de28fde3e500f4cefa8f48b97af48af.png" alt="">
|
||||
|
||||
上面这张图来自维基百科，显示的是26个字母在普通的英文文本中出现的概率。
|
||||
|
||||
只要我们能够获取足够长的密文进行分析的话，通过字母出现的频率，我们同样能够猜到相应的原始字母，这并不安全。所以，随着安全性需求的提高，人们有必要寻找一种容易将字母的自然频度隐蔽或均匀化，并使得统计分析足够安全可靠的加密方法。而希尔密码能基本满足这一要求，那么希尔密码是怎么做到这一点的呢？
|
||||
|
||||
## 希尔密码原理
|
||||
|
||||
我们先来看一下希尔密码的原理。根据百度百科的定义，希尔密码（Hill Cipher）是运用基本矩阵论原理的替换密码，由Lester S. Hill在1929年发明。每个字母当作26进制数字：A=0，B=1，C=2… ，把一串字母当成$n$维向量，和一个$n×n$的矩阵相乘，再将得出的结果和26进行模运算。
|
||||
|
||||
所以，希尔加密算法的基本思想是，通过线性变换将固定数量的明文字母转换为同样数量的密文字母，解密只要作一次逆变换就可以了，而密钥就是变换矩阵本身。
|
||||
|
||||
现在，我们再通过数学的方式来表达一下，希尔密码是如何通过三步来实现加密的。
|
||||
|
||||
第一步，设置加密矩阵$E$。
|
||||
|
||||
第二步，对照字母编码表（自行设定）得到数字，并把明文消息分割成大小为n的多个块：$v_{1},v_{2},…$，并且忽略空格。这里之所以忽略空格，是因为一般情况下密码传递的信息不会过于复杂。如果密码过于复杂，是可以分多次传递的。这里的n表示的密钥的阶数，密钥的阶数越高，也就是n越大的话，破译的难度也就越大，所需要的计算量也就越大。
|
||||
|
||||
第三步，每个消息块和加密矩阵$E$相乘：$Ev_{1}, Ev_{2},…$，并和26进行模运算，最后对照字母编码表得到密文。
|
||||
|
||||
同样，我们把这三步倒过来，就能实现解密了。
|
||||
|
||||
第一步，计算加密矩阵$E$的逆矩阵$D \equiv E^{-1}(\bmod 26)$。
|
||||
|
||||
第二步，对照字母编码表得到数字，把它和解密矩阵$D$相乘，并和26进行模运算。
|
||||
|
||||
第三步，对照编码表，得到原始明文。
|
||||
|
||||
这里你需要注意的是，加密矩阵很关键，它就是我们通常意义上所说的“密钥”，也就是打开密码的钥匙。
|
||||
|
||||
通过前面讲解的加密解密步骤，我们可以看出，希尔密码之所以很难被破译，是因为它设置了三道关卡：
|
||||
|
||||
1. 列矩阵的维度未知；
|
||||
1. 对应字母表的排列未知；
|
||||
1. 加密矩阵（或者说密钥）未知。
|
||||
|
||||
想要破解希尔密码，就需要同时获取到通过这三道关卡的钥匙，这谈何容易。
|
||||
|
||||
## 希尔密码实例
|
||||
|
||||
好了，原理都讲完了，现在我们通过一个例子来实际地看下希尔密码加密和解密的过程。
|
||||
|
||||
假设：A和B双方有一条重要消息要沟通，双方很早就建立了密钥沟通机制，每过一段时间都会更新密钥。在这次的密钥更新周期中，正确的密钥，也就是加密矩阵是一个3×3矩阵。
|
||||
|
||||
$$<br>
|
||||
E=\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
这一次A要给B的消息是“ILIKEBODYCOMBAT”，我们用之前的三步在A方先来加密：
|
||||
|
||||
第一步，定义加密矩阵，也就是刚才的$E$矩阵。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/9c/f1/9c7678f672c3fc46746b5504d336d0f1.png" alt="">
|
||||
|
||||
第二步，对照字母编码表得到数字：8、11、8、10、4、1、14、3、24、2、14、12、1、0、19。接下来，把明文消息分割成大小为3的5个块，也就是维度为3的5个列矩阵。
|
||||
|
||||
$$<br>
|
||||
v_{1}=\left[\begin{array}{c}<br>
|
||||
8 \\\<br>
|
||||
11 \\\<br>
|
||||
8<br>
|
||||
\end{array}\right], v_{2}=\left[\begin{array}{c}<br>
|
||||
10 \\\<br>
|
||||
4 \\\<br>
|
||||
1<br>
|
||||
\end{array}\right], v_{3}=\left[\begin{array}{c}<br>
|
||||
14 \\\<br>
|
||||
3 \\\<br>
|
||||
24<br>
|
||||
\end{array}\right], v_{4}=\left[\begin{array}{c}<br>
|
||||
2 \\\<br>
|
||||
14 \\\<br>
|
||||
12<br>
|
||||
\end{array}\right], v_{5}=\left[\begin{array}{c}<br>
|
||||
1 \\\<br>
|
||||
0 \\\<br>
|
||||
19<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第三步，将每个消息块和加密矩阵$E$相乘：
|
||||
|
||||
$$<br>
|
||||
E v_{1}=\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
8 \\\<br>
|
||||
11 \\\<br>
|
||||
8<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
320 \\\<br>
|
||||
360 \\\<br>
|
||||
467<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
8 \\\<br>
|
||||
22 \\\<br>
|
||||
25<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
E v_{2}=\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
10 \\\<br>
|
||||
4 \\\<br>
|
||||
1<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
157 \\\<br>
|
||||
204 \\\<br>
|
||||
283<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
1 \\\<br>
|
||||
22 \\\<br>
|
||||
23<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
E v_{3}=\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
14 \\\<br>
|
||||
3 \\\<br>
|
||||
24<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
180 \\\<br>
|
||||
470 \\\<br>
|
||||
691<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
24 \\\<br>
|
||||
2 \\\<br>
|
||||
15<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
E v_{4}=\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
2 \\\<br>
|
||||
14 \\\<br>
|
||||
12<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
360 \\\<br>
|
||||
370 \\\<br>
|
||||
458<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
22 \\\<br>
|
||||
6 \\\<br>
|
||||
16<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
E v_{5}=\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
1 \\\<br>
|
||||
0 \\\<br>
|
||||
19<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
25 \\\<br>
|
||||
203 \\\<br>
|
||||
305<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
25 \\\<br>
|
||||
21 \\\<br>
|
||||
19<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
最后，对照字母编码表得到密文：“IWZBWXBCGWGQZVT”。
|
||||
|
||||
B拿到这个密文后，使用三步来解密：
|
||||
|
||||
第一步，计算加密矩阵E的逆矩阵$D$：
|
||||
|
||||
$$<br>
|
||||
D \equiv\left[\begin{array}{ccc}<br>
|
||||
6 & 24 & 1 \\\<br>
|
||||
13 & 16 & 10 \\\<br>
|
||||
20 & 17 & 15<br>
|
||||
\end{array}\right]^{-1}(\bmod 26) \equiv\left[\begin{array}{ccc}<br>
|
||||
8 & 5 & 10 \\\<br>
|
||||
21 & 8 & 21 \\\<br>
|
||||
21 & 12 & 8<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第二步，对照字母编码表得到数字，把它和解密矩阵D相乘，并和26进行模运算，得到相应结果。
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{ccc}<br>
|
||||
8 & 5 & 10 \\\<br>
|
||||
21 & 8 & 21 \\\<br>
|
||||
21 & 12 & 8<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
8 \\\<br>
|
||||
22 \\\<br>
|
||||
25<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
424 \\\<br>
|
||||
869 \\\<br>
|
||||
632<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
8 \\\<br>
|
||||
11 \\\<br>
|
||||
8<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{ccc}<br>
|
||||
8 & 5 & 10 \\\<br>
|
||||
21 & 8 & 21 \\\<br>
|
||||
21 & 12 & 8<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
1 \\\<br>
|
||||
22 \\\<br>
|
||||
23<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
348 \\\<br>
|
||||
680 \\\<br>
|
||||
469<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
10 \\\<br>
|
||||
4 \\\<br>
|
||||
1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{ccc}<br>
|
||||
8 & 5 & 10 \\\<br>
|
||||
21 & 8 & 21 \\\<br>
|
||||
21 & 12 & 8<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
24 \\\<br>
|
||||
2 \\\<br>
|
||||
15<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
352 \\\<br>
|
||||
835 \\\<br>
|
||||
648<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
14 \\\<br>
|
||||
3 \\\<br>
|
||||
24<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{ccc}<br>
|
||||
8 & 5 & 10 \\\<br>
|
||||
21 & 8 & 21 \\\<br>
|
||||
21 & 12 & 8<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
22 \\\<br>
|
||||
6 \\\<br>
|
||||
16<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
366 \\\<br>
|
||||
846 \\\<br>
|
||||
662<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
2 \\\<br>
|
||||
14 \\\<br>
|
||||
12<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{ccc}<br>
|
||||
8 & 5 & 10 \\\<br>
|
||||
21 & 8 & 21 \\\<br>
|
||||
21 & 12 & 8<br>
|
||||
\end{array}\right]\left[\begin{array}{c}<br>
|
||||
25 \\\<br>
|
||||
21 \\\<br>
|
||||
19<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
495 \\\<br>
|
||||
1092 \\\<br>
|
||||
929<br>
|
||||
\end{array}\right] \bmod 26=\left[\begin{array}{c}<br>
|
||||
1 \\\<br>
|
||||
0 \\\<br>
|
||||
19<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
最后，B通过对照编码表，得到了原始明文：“ILIKEBODYCOMBAT”。
|
||||
|
||||
这里，你也许会问，密钥为什么用的是3×3的可逆矩阵？那是我为了例子方便而设置的，你完全可以设置更高阶的矩阵。就像之前说的，密钥的阶数越高，也就是$n$越大的话，破译的难度也就越大，所需要的计算量也就越大。
|
||||
|
||||
所以，从破译密码的角度来看，传统的密码有一个致命弱点，就是破译者可从统计出来的字符频率中找到规律，进而找出破译的突破口。尤其是在计算机技术高度发达的今天，破译的速度更快。而希尔密码算法则完全克服了这一缺陷，它通过采用线性代数中的矩阵乘法运算和逆运算，能够较好地抵抗频率分析，很难被攻破。
|
||||
|
||||
## 本节小结
|
||||
|
||||
这一节课的内容都和密码学有关，感觉像是搞谍战一样。但其实它的核心很简单，就是通过基础篇中学到的矩阵和逆矩阵的知识，来实现希尔密码。希尔密码的关键就是定义加密矩阵，或者说密钥、字母表排列方式和列矩阵的维度，通过线性变换将固定数量的明文字母转换为同样数量的密文字母，而解密则只要作一次逆变换就可以了。
|
||||
|
||||
当然，现实中还有更复杂的加密算法，其中最著名的，且用到线性代数的加密算法是AES，想必你平时也经常看到或用到过。AES是一个迭代的、对称密钥分组的密码，它可以使用128、192和256位密钥，并且用128、192和256位分组加密和解密数据，其中密钥长度与分组长度是独立的。
|
||||
|
||||
## 线性代数练习场
|
||||
|
||||
请你做一回“特工”，尝试使用希尔密码来给明文“MACHINELEARNING”做加密和解密。
|
||||
|
||||
>
|
||||
提醒：你可以自行定义加密矩阵、字母表排列方式和列矩阵的维度。加密矩阵可以使用之前介绍的3×3可逆矩阵，也可以使用其它n×n的可逆矩阵。
|
||||
|
||||
|
||||
欢迎在留言区晒出你的加密和解密过程，我会及时回复。同时，也欢迎你把这篇文章分享给你的朋友，一起讨论、学习。
|
||||
425
极客时间专栏/重学线性代数/应用篇/14 | 如何在深度学习中运用数值代数的迭代法做训练？.md
Normal file
425
极客时间专栏/重学线性代数/应用篇/14 | 如何在深度学习中运用数值代数的迭代法做训练？.md
Normal file
@@ -0,0 +1,425 @@
|
||||
<audio id="audio" title="14 | 如何在深度学习中运用数值代数的迭代法做训练？" controls="" preload="none"><source id="mp3" src="https://static001.geekbang.org/resource/audio/7b/cd/7b4f110f28feaa347ed455e7aa9a6ccd.mp3"></audio>
|
||||
|
||||
你好，我是朱维刚。欢迎你继续跟我学习线性代数，今天我要讲的内容是“数值线性代数的迭代法，以及如何在实践中运用迭代法求解线性方程组”。
|
||||
|
||||
大密度线性方程组的计算已经成为了世界上最快计算机的测试标准。2008年，IBM为美国能源部Los Alamos国家实验室建造了“Roadrunner”计算机系统，它的运算速度达到了1.026 petaflop/s（千万亿次/秒，petaflop是衡量计算机性能的一个重要单位，1 petaflop等于每秒钟进行1千万亿次的数学运算）。按摩尔定律计算，现在世界上最快的计算机已经达到了200 petaflop，我国也早就进入了世界前列，并有望实现1 exaflop/s（百亿亿次/秒），成为世界第一。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/25/f6/258e6b08673235b8a80ayy8da408a8f6.png" alt="">
|
||||
|
||||
可能你会有些疑惑，为什么我要在课程后期来讲数值线性代数呢？
|
||||
|
||||
那是因为数值线性代数是一门特殊的学科，是特别为计算机上进行线性代数计算服务的，可以说它是研究矩阵运算算法的学科，偏向算法实践与工程设计。有了之前基础知识的铺垫后，学习数值线性代数会更有效，而且它是可以直接运用在计算机科学中的，比如：在图像压缩中，使用奇异值分解（SVD）来节省内存；在深度学习中，使用共轭梯度来加速神经网络的收敛。
|
||||
|
||||
## 迭代方法说明
|
||||
|
||||
课程内容的前期一直都在用**直接法**来解线性方程组，比如高斯消元法。但在实践中，我们在面对复杂场景时，更多的会使用**迭代法**来求解（也就是所谓的间接法），因为很多场景会用到大型稀疏矩阵。所以，我打算在这里讲讲机器学习中的迭代法应用。这里需要注意，不是说直接法不重要，直接法解决了许多相对简单的问题，也是其他方法的基础。
|
||||
|
||||
现在我就来说一说什么是迭代法？
|
||||
|
||||
我们还是通过线性方程组$Ax=b$来看看。在这里我们分解$A$，使得$A=S-T$，代入等式后得出：$Sx=Tx+b$（等式①）。
|
||||
|
||||
按这样的方式持续下去，通过迭代的方式来解$Sx$。这就类似于把复杂问题层层分解和简化，最终使得这个迭代等式成立：$Sx_{k+1}=Tx_{k}+b$（等式②）。
|
||||
|
||||
更具体一点来说，我们其实是从$x_{0}$开始，解$Sx_{1}=Tx_{0}+b$。然后，继续解$Sx_{2}=Tx_{1}+b$。一直到$x_{k+1}$非常接近$x_{k}$时，又或者说残余值$r_{k}=b-Ax_{k}$接近$0$时，迭代停止。由于线性方程组的复杂程度不同，这个过程经历几百次的迭代都是有可能的。所以，迭代法的目标就是**比消元法更快速地逼近真实解**。
|
||||
|
||||
那么究竟应该如何快速地逼近真实解呢？
|
||||
|
||||
这里，$A=S-T$，$A$的分解成了关键，也就是说$A$的分解目标是**每步的运算速度和收敛速度都要快**。每步的运算速度取决于$S$，而收敛速度取决于“错误”(error)，$e_{k}$，这里的错误 $e_{k}$是$x-x_{k}$，也就是说$x$和$x_{k}$的差应该快速逼近0，我们把错误表示成这样：$e_{k+1}=S^{-1}Te_{k}$（等式③）。
|
||||
|
||||
它是等式②和①差后得出的结果，迭代的每一步里，错误都会被$S^{-1}T$乘，如果$S^{-1}T$越小，那逼近0的速度就更快。在极端分解情况下，$S=A$，$T=0$，那$Ax=b$又回来了，第一次迭代就能完成收敛，其中$S^{-1}T$等于0。
|
||||
|
||||
但是，这一次迭代的成本太高，我们回到了非迭代方式的原点。所以，你也知道，鱼和熊掌不能兼得，$S$的选择成为了关键。那我们要如何在每一次迭代的速度和快速收敛之间做出平衡呢？我给你$S$选择的几种常见方法：
|
||||
|
||||
1. 雅可比方法（Jacobi method）：$S$取$A$的对角部分。
|
||||
1. 高斯-赛德尔方法（Gauss-Seidel）：$S$取$A$的下三角部分，包含对角。
|
||||
1. ILU方法（Incomplete LU）：$S=L$估计乘$U$估计。
|
||||
|
||||
## 雅可比方法实践
|
||||
|
||||
总体介绍了迭代法理论之后，我们就进入迭代法运用的实践环节。
|
||||
|
||||
首先，我们先来试试使用雅可比方法解线性方程组，雅克比迭代法是众多迭代法中比较早且较简单的一种。所以，作为迭代法的实践开篇比较合适。让我们设一个2×2的线性方程组：
|
||||
|
||||
$$<br>
|
||||
Ax=b<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
\left\{\begin{array}{c}<br>
|
||||
2 u-v=4 \\\<br>
|
||||
-u+2 v=-2<br>
|
||||
\end{array}\right.<br>
|
||||
$$
|
||||
|
||||
我们很容易就能得出这个方程组的解如下。
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u \\\<br>
|
||||
v<br>
|
||||
\end{array}\right]=\left[\begin{array}{l}<br>
|
||||
2 \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
现在我们就用雅可比方法来看看怎么解这个方程组：
|
||||
|
||||
首先，我们把线性方程组转换成矩阵形式。
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{cc}<br>
|
||||
2 & -1 \\\<br>
|
||||
-1 & 2<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
4 \\\<br>
|
||||
-2<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
接着，把A的对角线放在等式左边，得出$S$矩阵。
|
||||
|
||||
$$<br>
|
||||
S=\left[\begin{array}{ll}<br>
|
||||
2 & 0 \\\<br>
|
||||
0 & 2<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
其余部分移到等式右边，得出$T$矩阵。
|
||||
|
||||
$$<br>
|
||||
T=\left[\begin{array}{ll}<br>
|
||||
0 & 1 \\\<br>
|
||||
1 & 0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
于是，雅可比迭代就可以表示成下面这样的形式。
|
||||
|
||||
$$<br>
|
||||
\mathrm{S} x_{k+1}=T x_{k}+\mathrm{b}<br>
|
||||
$$
|
||||
|
||||
$$<br>
|
||||
\left\{\begin{array}{l}<br>
|
||||
2 u_{k+1}=v_{k}+4 \\\<br>
|
||||
2 v_{k+1}=u_{k}-2<br>
|
||||
\end{array}\right.<br>
|
||||
$$
|
||||
|
||||
现在是时候进行迭代了，我们从$u_{0}=v_{0}=0$开始。
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{0} \\\<br>
|
||||
v_{0}<br>
|
||||
\end{array}\right]=\left[\begin{array}{l}<br>
|
||||
0 \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第一次迭代后，我们得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{1} \\\<br>
|
||||
v_{1}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
2 \\\<br>
|
||||
-1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第二次迭代后得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{2} \\\<br>
|
||||
v_{2}<br>
|
||||
\end{array}\right]=\left[\begin{array}{l}<br>
|
||||
\frac{3}{2} \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第三次迭代后得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{3} \\\<br>
|
||||
v_{3}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
2 \\\<br>
|
||||
-\frac{1}{4}<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第四次迭代后得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{4} \\\<br>
|
||||
v_{4}<br>
|
||||
\end{array}\right]=\left[\begin{array}{l}<br>
|
||||
\frac{15}{8} \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第五次迭代后，我们得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{5} \\\<br>
|
||||
v_{5}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
2 \\\<br>
|
||||
-\frac{1}{16}<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
经过五次迭代后发现收敛，因为它的结果接近真实解。
|
||||
|
||||
$$<br>
|
||||
真实解<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
2 \\\<br>
|
||||
0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
现在，再来看一下错误等式，$\mathrm{Se}**{k+1}=T e**{k}$，我们把$S$和$T$代入等式，得出：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{ll}<br>
|
||||
2 & 0 \\\<br>
|
||||
0 & 2<br>
|
||||
\end{array}\right] e_{k+1}=\left[\begin{array}{ll}<br>
|
||||
0 & 1 \\\<br>
|
||||
1 & 0<br>
|
||||
\end{array}\right] e_{k}<br>
|
||||
$$
|
||||
|
||||
计算$S$的逆矩阵和$T$相乘$S^{-1}T$得出：
|
||||
|
||||
$$<br>
|
||||
e_{k+1}=\left[\begin{array}{cc}<br>
|
||||
0 & \frac{1}{2} \\\<br>
|
||||
\frac{1}{2} & 0<br>
|
||||
\end{array}\right] e_{k}<br>
|
||||
$$
|
||||
|
||||
这里，$S$的逆矩阵和T相乘$S^{-1}T$有特征值$\frac{1}{2}$和$-\frac{1}{2}$，所以，它的谱半径是$\rho(B)=\frac{1}{2}$。这里的**谱半径**是用来控制收敛的，所以非常重要。谱半径从数学定义上是：矩阵（或者有界线性算子的谱半径）是指其特征值绝对值集合的上确界。这个概念是不是很难理解？具体谱半径的概念你可以查互联网来获取，为了方便你理解，这里我还是用数学方法来简单表达一下。
|
||||
|
||||
$$<br>
|
||||
B=S^{-1} T=\left[\begin{array}{ll}<br>
|
||||
0 & \frac{1}{2} \\\<br>
|
||||
\frac{1}{2} & 0<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
通过$S$的逆矩阵和$T$相乘$S^{-1}T$，我们得到：$|\lambda|_{\max }=\frac{1}{2}$，以及：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{cc}<br>
|
||||
0 & \frac{1}{2} \\\<br>
|
||||
\frac{1}{2} & 0<br>
|
||||
\end{array}\right]^{2}=\left[\begin{array}{cc}<br>
|
||||
\frac{1}{4} & 0 \\\<br>
|
||||
0 & \frac{1}{4}<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
这里的特征值$\frac{1}{2}$非常小，所以10次迭代后，错误就很低了，即$\frac{1^{10}}{2}=\frac{1}{1024}$。而如果特征值是0.99或者0.999，那很显然迭代次数就要多得多，也就是说需要更多时间来做运算。
|
||||
|
||||
## 高斯-赛德尔方法实践
|
||||
|
||||
现在我们再来看下高斯-赛德尔方法，高斯-赛德尔迭代可以**节约存储**和**加速迭代**，每迭代一次只需一组存储单元，而雅可比迭代需要两组单元。
|
||||
|
||||
$S$取$A$的下三角部分，还是使用之前雅可比方法中的例子，我们得出方程组：
|
||||
|
||||
$$<br>
|
||||
\left\{\begin{array}{c}<br>
|
||||
u_{k+1}=\frac{1}{2} v_{k}+2 \\\<br>
|
||||
v_{k+1}=\frac{1}{2} u_{k+1}-1<br>
|
||||
\end{array}\right.<br>
|
||||
$$
|
||||
|
||||
这里有一个比较大的变化，那就是$u_{k}$消失了，通过$v_{k}$，我们可以直接得到$u_{k+1}$和$v_{k+1}$，这样有什么好处呢？两大好处是显而易见的，就是**节约存储**和**加速迭代**。
|
||||
|
||||
接下来，我们从$u_{0}=0$，$v_{0}=-1$来测试一下迭代。
|
||||
|
||||
第一次迭代后，我们得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{1} \\\<br>
|
||||
v_{1}<br>
|
||||
\end{array}\right]=\left[\begin{array}{c}<br>
|
||||
\frac{3}{2} \\\<br>
|
||||
\frac{-1}{4}<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第二次迭代后得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{2} \\\<br>
|
||||
v_{2}<br>
|
||||
\end{array}\right]=\left[\begin{array}{l}<br>
|
||||
\frac{15}{8} \\\<br>
|
||||
\frac{-1}{16}<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
第三次迭代后得到：
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{l}<br>
|
||||
u_{3} \\\<br>
|
||||
v_{3}<br>
|
||||
\end{array}\right]=\left[\begin{array}{r}<br>
|
||||
\frac{63}{32} \\\<br>
|
||||
-\frac{1}{64}<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
经过三次迭代后发现收敛，因为第三次迭代后的结果接近真实解。
|
||||
|
||||
错误经过计算分别是$-1, \frac{-1}{4}, \frac{-1}{16}, \frac{-1}{64}$，和刚才使用雅可比方法得出的错误$2, \frac{1}{2}, \frac{1}{8}, \frac{1}{32}$。比较后我们可以发现，无论是迭代次数还是收敛速度方面，高斯-赛德尔方法比雅可比方法速度快、精确度也高得多。
|
||||
|
||||
## 逐次超松弛方法
|
||||
|
||||
最后，我们在高斯-赛德尔方法上做个小调整，在迭代中引入一个参数“omega”，$ω$，即超松弛因子。然后选择一个合适的$ω$，使得$S^{-1}T$的谱半径尽可能小，这个方法就叫做逐次超松弛方法（Successive over-relaxation method，简称SOR）。
|
||||
|
||||
SOR方法的方程是：$ωAx=ωb$，矩阵$S$有$A$的对角线，对角线下是$ωA$，等式右边$T$是$S-ωA$，于是，我们还是使用之前雅可比方法中的例子，得到SOR方程组如下。
|
||||
|
||||
$$\left\{\begin{array}{c}<br>
|
||||
2 u_{k+1}=(2-2 \omega) u_{k}+\omega v_{k}+4 \omega \\\<br>
|
||||
-\omega u_{k+1}+2 v_{k+1}=(2-2 \omega) v_{k}-2 \omega<br>
|
||||
\end{array}\right.$$
|
||||
|
||||
是不是看起来更复杂了？
|
||||
|
||||
没关系，其实它只是在我们眼中看起来复杂，对计算机来说是没区别的。对SOR来说，只是多了一个$ω$，而$ω$选择越好就越快。具体$ω$的选择，以及迭代的过程就不赘述了，我给你一个小提示，你可以在“$ω$大于1”和“$ω$小于1”两种情况下来多选择几个$ω$进行尝试，最后你应该会得到结论：
|
||||
|
||||
1. 在$ω$大于1时，$ω$越大，迭代的次数就越多，收敛速度就越慢，$ω$接近1时，迭代的次数越小，收敛速度越快。
|
||||
1. 在$ω$小于1时，$ω$越小，迭代的次数就越多，收敛速度就越慢，$ω$接近1时，迭代的次数越小，收敛速度越快。
|
||||
|
||||
所以，SOR迭代法的关键就是$ω$的选择，它可以被看作是高斯-赛德尔法的扩充。
|
||||
|
||||
雅可比法、高斯-赛德尔法，以及SOR迭代法都是定常迭代法。接下来我讲一下和定常迭代法不同的另一类方法，也是实践中用的比较多的方法——**共轭梯度法**（Conjugate gradient），它属于Krylov子空间方法。简单来说，Krylov子空间方法是一种 “降维打击” 手段，是一种牺牲精度换取速度的方法。
|
||||
|
||||
## 共轭梯度法
|
||||
|
||||
要讲共轭梯度法，我们要先解释一下“共轭”，共轭就是按一定的规律相配的一对，通俗点说就是孪生。“轭”是牛拉车用的木头，那什么是共轭关系呢？同时拉一辆车的两头牛，就是共轭关系。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/10/c1/10a99c815d4961ed00c8f4a8693c65c1.jpg" alt="">
|
||||
|
||||
我们根据这个定义再来解释一下共轭方向，向量$p, q \in R$， 若满足条件$pAq=0$， 则称$p$和$q$关于$A$是共轭方向，或者$p$和$q$关于$A$共轭。有了共轭和共轭方向的概念后，再来看共轭梯度法就简单多了。共轭梯度法的出现不仅是为了解决梯度下降法的收敛速度慢，而且也避免了牛顿法需要存储和计算黑塞矩阵（Hessian Matrix）并求逆的缺点。
|
||||
|
||||
现在来看看共轭梯度算法，设$Ax=b$，其中$A$是一个实对称正定矩阵。
|
||||
|
||||
首先，我们设初始值$x_{0}$为$0$或者一个估计值，来计算$r_{0}:=b-A x_{0}$。如果$r_{0}$非常小，那$x_{0}$就是结果，如果不是就继续。
|
||||
|
||||
接下来设$p_{0}:=r_{0}$，$k:=0$。现在我们开始迭代循环。
|
||||
|
||||
a.计算$\alpha_{k}$ 。
|
||||
|
||||
$$\alpha_{k}:=\frac{r_{k}^{T} r_{k}}{p_{k}^{T} A p_{k}}$$
|
||||
|
||||
b.计算$x_{k+1}$ 。
|
||||
|
||||
$$x_{k+1}:=x_{k}+\alpha_{k} p_{k}$$
|
||||
|
||||
c.计算$r_{k+1}$。
|
||||
|
||||
$$r_{k+1}:=r_{k}-\alpha_{k} A p_{k}$$
|
||||
|
||||
d.如果$r_{k+1}$非常小，循环结束，如果不是就继续。
|
||||
|
||||
e.计算$β_{k}$ 。
|
||||
|
||||
$$\beta_{k}:=\frac{r_{k+1}^{T} r_{k+1}}{r_{k}^{T} r_{k}}$$
|
||||
|
||||
f.计算$p_{k+1}$ 。
|
||||
|
||||
$$p_{k+1}:=r_{k+1}+\beta_{k} p_{k}$$
|
||||
|
||||
g.$k:=k+1$。
|
||||
|
||||
1. 返回结果$x_{k+1}$。
|
||||
|
||||
从算法中我们可以看出，共轭梯度法的优点是**存储量小**和**具有步收敛性**。如果你熟悉MATLAB，就会发现共轭梯度法的实现超级简单，只需要短短十几行代码（下方代码来自于MATLAB/GNU Octave的例子）。
|
||||
|
||||
```
|
||||
function x = conjgrad(A, b, x)
|
||||
r = b - A * x;
|
||||
p = r;
|
||||
rsold = r' * r;
|
||||
|
||||
|
||||
for i = 1:length(b)
|
||||
Ap = A * p;
|
||||
alpha = rsold / (p' * Ap);
|
||||
x = x + alpha * p;
|
||||
r = r - alpha * Ap;
|
||||
rsnew = r' * r;
|
||||
if sqrt(rsnew) < 1e-10
|
||||
break;
|
||||
end
|
||||
p = r + (rsnew / rsold) * p;
|
||||
rsold = rsnew;
|
||||
end
|
||||
end
|
||||
|
||||
```
|
||||
|
||||
## 机器学习中的共轭梯度
|
||||
|
||||
共轭梯度法经常被用在训练神经网络中，在实践中已经证明，它是比**梯度下降**更有效的方法，因为就像刚才讲的，它不需要计算黑塞矩阵。那我现在就来讲一讲，使用共轭梯度法的神经网络训练过程。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/b4/a1/b45e4b977f7c21c747844c383985b1a1.png" alt="">
|
||||
|
||||
在整个训练过程中，**参数改进**是重点，当然这也是所有神经网络训练的重点。这个过程是通过计算共轭梯度的训练方向，然后计算训练速率来实现的。在共轭梯度训练算法中，搜索是按共轭方向进行的，也就是说，训练方向是共轭的。所以，收敛速度比梯度下降要快。
|
||||
|
||||
现在我们来看训练方向的计算方法。首先，我们设置训练方向向量为$d$，然后，定义一个初始参数向量$w^{0}$，以及一个初始训练方向向量$d^{0}=-g^{0}$，于是，共轭梯度法构造出的训练方向可以表示成：$d^{i+1}=g^{i+1}+d^{i} \cdot \gamma^{i}$。
|
||||
|
||||
其中，$g$是梯度向量，$γ$是共轭参数。参数通过这个表达式来更新和优化。通常训练速率$\eta$可使用单变量函数优化方法求得。
|
||||
|
||||
$$w^{i+1}=w^{i}+d^{i} \cdot \eta^{i}$$
|
||||
|
||||
## 本节小结
|
||||
|
||||
好了，到这里数值线性代数的迭代法这一讲就结束了，最后我再总结一下前面讲解的内容。
|
||||
|
||||
首先，我先解释了数值线性代数，接着再整体讲解了迭代方法。然后，举了一个线性方程组的例子，运用迭代法中的几个比较著名的实践方法：雅可比方法、高斯-赛德尔方法，以及逐次超松弛方法，来解这个线性方程组。最后，我把共轭梯度法用在了深度学习的神经网络训练中。
|
||||
|
||||
希望你能在了解了数值线性代数，以及迭代法后，更多地在计算机科学领域中，运用迭代法做矩阵运算。如果有兴趣，你也可以学习其它在实践中使用的迭代法。
|
||||
|
||||
## 线性代数练习场
|
||||
|
||||
练习时刻到了，这次继续使用第一篇线性方程组里的例子，你可以挑选任意一个迭代法来求解这个线性方程组。
|
||||
|
||||
假设，一个旅游团由孩子和大人组成，去程时他们一起坐大巴，每个孩子的票价3元，大人票价3.2元，总共花费118.4元。回程时一起做火车，每个孩子的票价3.5元，大人票价3.6元，总共花费135.2元。请问这个旅游团中有多少孩子和大人？
|
||||
|
||||
设小孩人数为$x_{1}$，大人人数为$x_{2}$，于是我们得到了一个方程组：
|
||||
|
||||
$$\left\{\begin{array}{c}<br>
|
||||
3 x_{1}+3.2 x_{2}=118.4 \\\<br>
|
||||
3.5 x_{1}+3.6 x_{2}=135.2<br>
|
||||
\end{array}\right.$$
|
||||
|
||||
这个方程组的解是：
|
||||
|
||||
$$\left\{\begin{array}{l}<br>
|
||||
x_{1}=16 \\\<br>
|
||||
x_{2}=22<br>
|
||||
\end{array}\right.$$
|
||||
|
||||
你可以计算一下多少次迭代后它能收敛，也就是逼近真实解？以及它的错误$e$又分别是多少？
|
||||
|
||||
欢迎在留言区晒出的你的运算过程和结果。如果有收获，也欢迎你把这篇文章分享给你的朋友。
|
||||
137
极客时间专栏/重学线性代数/应用篇/15 | 如何从计算机的角度来理解线性代数？.md
Normal file
137
极客时间专栏/重学线性代数/应用篇/15 | 如何从计算机的角度来理解线性代数？.md
Normal file
@@ -0,0 +1,137 @@
|
||||
<audio id="audio" title="15 | 如何从计算机的角度来理解线性代数？" controls="" preload="none"><source id="mp3" src="https://static001.geekbang.org/resource/audio/5a/e3/5afe62a286708b19cfde3702da0064e3.mp3"></audio>
|
||||
|
||||
你好，我是朱维刚。欢迎你继续跟我学习线性代数，今天我要讲的内容是“如何从计算机的角度来理解线性代数”。
|
||||
|
||||
基础和应用篇整体走了一圈后，最终我们还是要回归到一个话题——从计算机的角度来理解线性代数。或者更确切地说，如何让计算机在保证计算精度和内存可控的情况下，快速处理矩阵运算。在数据科学中，大部分内容都和矩阵运算有关，因为几乎所有的数据类型都能被表达成矩阵，比如：结构化数据、时序数据、在Excel里表达的数据、SQL数据库、图像、信号、语言等等。
|
||||
|
||||
线性代数一旦和计算机结合起来，需要考虑的事情就多了。你还记得开篇词中我讲到的四个层次的最后一层——“能够踏入大规模矩阵计算的世界”吗？当我们面对大规模矩阵的时候，计算机的硬件指标就需要考虑在内了，这也是硬性的限制条件。在碰到大规模矩阵的时候，这些限制条件会被放大，所以精度、内存、速度和扩展这四点是需要你思考的。
|
||||
|
||||
1. 精度：计算机的数字计算是有有限精度的，这个想必你能理解，当遇到迭代计算的情况下，四舍五入会放大很小的误差；
|
||||
1. 内存：一些特殊结构的矩阵，比如包含很多0元素的矩阵，可以考虑优化内存存储方式；
|
||||
1. 速度：不同的算法、并行执行、以及内存数据移动的耗时，这些都和速度有关；
|
||||
1. 扩展：当单机内存不够时，你在考虑横向扩展的同时，还要考虑如何分片，也就是如何分布矩阵运算的算力。
|
||||
|
||||
接下来，我就从这几点深入讲解一下。
|
||||
|
||||
## 精度
|
||||
|
||||
首先是精度，我们先从计算机如何存储数字的角度入手，来做一个练习。你可以执行一下这个Python代码，想想会发生什么情况呢？
|
||||
|
||||
```
|
||||
def f(x):
|
||||
if x <= 1/2:
|
||||
return 2 * x
|
||||
if x > 1/2:
|
||||
return 2*x - 1
|
||||
|
||||
x = 1/10
|
||||
|
||||
for i in range(80):
|
||||
print(x)
|
||||
x = f(x)
|
||||
|
||||
```
|
||||
|
||||
这个结果可能会和你想象的有很大的不同。
|
||||
|
||||
其实数学和计算机之间存在很大的不同，数学是连续的、无穷的，而计算机是离散的、有限的。就像这个练习，计算机存储的数字精度是有限的，而很小的误差通过很多次的迭代会累加，最终放大成比较大的错误。一个惨痛的例子就是1996年欧洲航天局的Ariane 5号火箭的发射失败，最后发现问题出在操作数误差上，也就是64位数字适配16位空间发生的整数溢出错误。
|
||||
|
||||
那我们该怎么来理解**数学是连续的，而计算机是离散的**呢？举个简单的例子，我们来看下数学中的区间表达[1,2]，这个形式就是连续的；但如果在计算机中以双精度来表达同样的东西，则会是这样的离散形式：
|
||||
|
||||
$$<br>
|
||||
1,1+2^{-52}, 1+2 \times 2^{-52}, 1+3 \times 2^{-52}, \ldots, 2<br>
|
||||
$$
|
||||
|
||||
于是，我们就引出了一个计算机领域精度计算的概念——机械最小值（EPSILON），对于双精度来说，IEEE标准指定机械最小值是：$ε=2^{-53}$。
|
||||
|
||||
## 内存
|
||||
|
||||
刚才我们看的是数字精度，现在我们接着来看看矩阵的存储方式。我们都知道，内存是有限的，所以，当你面对大矩阵时，千万不要想着把矩阵所有的元素都存储起来。解决这个问题的一个最好方式就是只存储非零元素，这种方式就叫做稀疏存储。稀疏存储和稀疏矩阵是完美匹配的，因为稀疏矩阵大部分的元素都是零。
|
||||
|
||||
我们来看一个机器学习中比较简单的稀疏矩阵的例子：Word Embedding中的One-Hot编码。One-Hot编码就是给句子中的每个字分别用一个0或1编码，一个句子中有多少个字，就有多少维度。这样构造出来的矩阵是很大的，而且是稀疏矩阵。比如：“重学线性代数”这六个字，通过One-Hot编码，就能表达成下面这样的形式。
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{llllll}<br>
|
||||
1 & 0 & 0 & 0 & 0 & 0 \\\<br>
|
||||
0 & 1 & 0 & 0 & 0 & 0 \\\<br>
|
||||
0 & 0 & 1 & 0 & 0 & 0 \\\<br>
|
||||
0 & 0 & 0 & 1 & 0 & 0 \\\<br>
|
||||
0 & 0 & 0 & 0 & 1 & 0 \\\<br>
|
||||
0 & 0 & 0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
所以，一般稀疏矩阵的大致形态如下图所示。
|
||||
|
||||
$$<br>
|
||||
\left[\begin{array}{llllllll}<br>
|
||||
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\<br>
|
||||
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\<br>
|
||||
0 & 0 & 2 & 0 & \frac{1}{2} & 0 & 0 & 0 \\\<br>
|
||||
0 & 0 & 0 & 1 & 0 & 0 & 2 & 0 \\\<br>
|
||||
0 & 0 & \frac{1}{2} & 0 & 4 & 0 & 0 & 0 \\\<br>
|
||||
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\\<br>
|
||||
0 & 0 & 0 & 1 & 2 & 0 & 1 & 0 \\\<br>
|
||||
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1<br>
|
||||
\end{array}\right]<br>
|
||||
$$
|
||||
|
||||
也有特殊类型的结构化矩阵，比如：对角线矩阵、三对角矩阵、海森堡矩阵等等，它们都有自己的特殊稀疏表达，也都通常被用来减少内存存储和计算量。可以说，数值线性代数在运算方面更多地聚焦在稀疏性上。
|
||||
|
||||
## 速度
|
||||
|
||||
接下来我们再来看速度。速度涉及到很多方面，比如计算复杂度、单指令多数据矢量运算、存储类型和网络等。
|
||||
|
||||
### 计算复杂度
|
||||
|
||||
算法通常由计算复杂度表达，同一问题可用不同算法解决，而一个算法的质量优劣将影响到算法乃至程序的效率。算法分析的目的在于选择合适的算法或者优化算法。
|
||||
|
||||
那么我们该如何评价一个算法呢？主要就是从时间复杂度和空间复杂度来综合考虑的。从矩阵计算的角度看，计算复杂度的考虑是很有必要的。简单的计算，我们可以用直接法来计算，而有的比较复杂的计算就要用间接迭代法了。特别是在很多场景中会用到的大型稀疏矩阵，这些计算复杂度都是不同的。我在第4篇“[解线性方程组](https://time.geekbang.org/column/article/269448)”中提到了迭代法，你可以回顾一下，同时你也可以参考[上一节课](https://time.geekbang.org/column/article/279476)的迭代法应用细节。
|
||||
|
||||
### 单指令多数据矢量运算
|
||||
|
||||
现代CPU和GPU都支持在同一时间以同步方式执行同一条指令。
|
||||
|
||||
举个例子来说，一个向量中的4个浮点数的指数幂运算可以同时执行，从而极大地提高运算效率，这类单指令多数据矢量运算处理就叫做SIMD（Single Instruction Multiple Data）。这在矩阵运算中非常重要，虽然我们不用太关心底层的实现，但如果你可以了解一些矩阵运算的包和库，那么你就可以在实际的开发中直接使用它们了，其中比较出名的就是python的NumPy，以及BLAS（Basic Linear Algebra Subprograms）和LAPACK。
|
||||
|
||||
>
|
||||
LAPACK是用Fortran写的，这里也纪念一下Fortran创始人约翰·巴库斯(John W. Backus)，不久前在美国俄勒冈州的家中去世，享年82岁。
|
||||
|
||||
|
||||
### 存储类型和网络
|
||||
|
||||
计算机存储类型有很多，比如：缓存、内存、机械盘、SSD，这些不同存储媒介的存储延迟都是不一样的；在网络方面，不同IDC、地区、地域的传输速率也是不同的。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/c9/30/c9a520f0c2933951be9240cd3e36ee30.png" alt="">
|
||||
|
||||
上图中的这些数据是所有程序员必须要知道的，因为毕竟各类计算机资源是有限的，在做解决方案时，必须综合考虑存储和网络的性能和成本来做组合，最终达到最大的产出投入比。这些数据来自[GitHub](https://colin-scott.github.io/personal_website/research/interactive_latency.html)，可以调整年限值来观察每年的动态变换。
|
||||
|
||||
## 扩展
|
||||
|
||||
最后我再来讲一下扩展，扩展分为单机多核、多处理器的垂直方向扩展，以及多节点平行方向扩展。
|
||||
|
||||
当我们想要利用多处理器能力来处理大型矩阵运算的时候，传统的方法就是把大矩阵分解成小矩阵块。比如：一台拥有4个处理器的服务器，现在有这样的两个6×6矩阵$A$和$B$要做乘运算。
|
||||
|
||||
<img src="https://static001.geekbang.org/resource/image/32/bb/32cec179f0604fe3f41ca87b9aa324bb.png" alt=""><img src="https://static001.geekbang.org/resource/image/72/48/72d70b50b3d110f0e404fa29be400548.png" alt="">
|
||||
|
||||
我们可以把这两个矩阵分别分成四块，每块都是3×3矩阵，例如：$A11、A12、A21$和$A22、B11、B12、B21$和$B22$，$A$和$B$的乘运算就是把各自分好的块分配到各处理器上做并行处理，处理后的结果再做合并。
|
||||
|
||||
比如：处理器P1处理$A11$和$B11$，处理器P2处理$A12$和$B12$，处理器P3处理$A21$和$B21$，处理器P4处理$A22$和$B22$。于是，4个处理器分别处理$A$和$B$的乘计算来获取结果：$C11、C12、C21$和$C22$。拿$C11$来看，$C11=A11×B11+A12×B21$，虽然$B21$不在P2中，但可以从P3传递过来再计算。
|
||||
|
||||
当计算机垂直扩展到极限后，就需要考虑扩展到多节点计算了，其实原理也是一样的，不一样的是需要从应用层面来设计调度器，来调度不同的计算机节点来做计算。
|
||||
|
||||
## 本节小结
|
||||
|
||||
线性代数运用在计算机科学中就是数值线性代数，它是一门特殊的学科，是特别为计算机上进行线性代数计算服务的，可以说它是研究矩阵运算算法的学科。这里，你需要掌握很重要的一个点就是：数学和计算机之间存在很大的不同，数学是连续的、无穷的，而计算机是离散的、有限的。
|
||||
|
||||
所以，从计算机角度来执行矩阵运算，需要考虑很多方面：精度、内存、速度和扩展，这样，你在做解决方案时，又或者在写程序时，才能在计算机资源有限的情况下做到方案或程序的最优化，也可以避免类似1996年欧洲航天局的Ariane 5号火箭发射失败的这类错误。
|
||||
|
||||
## 线性代数练习场
|
||||
|
||||
我想让最后一篇的练习成为一个知识点的补充。
|
||||
|
||||
针对矩阵高性能并行计算，我前面在“扩展”一模块讲的是一般传统方法，也就是把大矩阵分解成小矩阵块，利用多处理器能力来处理大型矩阵运算。现在，请你研究一下如何用Cannon算法解决这个问题？
|
||||
|
||||
Cannon算法是一种存储效率很高的算法，也是对传统算法的改进，目标就是减少分块矩阵乘法的存储量。而且你也可以把它看成是MPI编程的一个例子。
|
||||
|
||||
欢迎在留言区分享你的研究成果，大家一起探讨。同时，也欢迎你把这篇文章分享给你的朋友，一起讨论、学习。
|
||||
11
极客时间专栏/重学线性代数/应用篇/强化通关 | 线性代数水平测试20题.md
Normal file
11
极客时间专栏/重学线性代数/应用篇/强化通关 | 线性代数水平测试20题.md
Normal file
@@ -0,0 +1,11 @@
|
||||
<audio id="audio" title="强化通关 | 线性代数水平测试20题" controls="" preload="none"><source id="mp3" src="https://static001.geekbang.org/resource/audio/8f/29/8f371dd588b2a4a8d98fb72c70a49b29.mp3"></audio>
|
||||
|
||||
你好，我是朱维刚。
|
||||
|
||||
《重学线性代数》这门课程就正式完结了，很感谢你一直以来的认真学习和支持！
|
||||
|
||||
为了帮助你检验自己的学习效果，我特别给你准备了一套结课测试题。这套测试题一共有 20 道题目，都是单选题，其中大部分都来自历届的自学考和考研的线性代数题，来重新回顾一下当年吧。试卷满分 100 分，我建议你在 45 分钟以内完成。
|
||||
|
||||
还等什么，快点击下面的按钮开始测试吧！
|
||||
|
||||
[<img src="https://static001.geekbang.org/resource/image/28/a4/28d1be62669b4f3cc01c36466bf811a4.png" alt="">](http://time.geekbang.org/quiz/intro?act_id=208&exam_id=583)
|
||||
Reference in New Issue
Block a user