https://web.archive.org/web/20131017050335/http://blog.renren.com/share/277469879/16492231759
9301 三维空间中二维网格的迭代结果不是环面来源： 王赟 ^_^Maigo 的日志 昨天， @ 张文琦(228254913) 同学的一篇日志又把椭圆迭代问题推向了高潮。这篇日志把平面中的多边形推广到了三维空间中的二维网格。作者得到的结论是，二维网格的迭代结果是环面（ torus ）。同时作者注意到，有时候环面上会打一个「结」，这个结需要大量的迭代才能解开。
我对这个「结」感到非常好奇，于是继续用上一篇日志中的傅里叶分析方法对二维网格的迭代进行了分析。结果发现：这样的「结」其实有两个；它们并不只是「有时」会出现，而是必然会出现；而且，无论迭代多少次，这两个「结」都不会消失。二维网格的迭代结果并不是环面，而是一种叫作 elliptic translational celestial 的曲面，它与环面最显著的区别就是有两个「结」。
photo-media
一、三维空间中二维网格迭代问题的描述
与上一篇日志一样，我用小括号表示迭代次数，中括号表示点的下标。
首先定义二维网格。在三维空间中取 n*m 个点 P[i][j] = (x[i][j], y[i][j], z[i][j])，其中 0<=i<n ， 0<=j<m。另设 P[n][j] = P[0][j] ， P[i][m] = P[i][0]。将所有一个下标相同、另一个下标相差 1 的点对相连，就组成了一个网格。这个网格的每一个「孔」是一个（空间）四边形，四边形 [i][j] 的四个顶点为 P[i][j]、P[i+1][j]、P[i+1][j+1]、P[i][j+1]。
迭代过程如下：首先取定初始网格的各个顶点 P(0)[i][j] ，并平移使得各个顶点的平均值位于原点。每次迭代时，取各个四边形的中心为新网格的顶点，即 P(k)[i][j] = ( P[i][j] + P[i+1][j] + P[i+1][j+1] + P[i][j+1] ) / 4 ，然后缩放使得各点的 x、y、z 坐标的平方和均为 1。
问：二维网格经过上述迭代过程，收敛成什么形状？
二、对二维网格形状的朴素思考
读到这里，你脑中想象出来的二维网格也许是分布在一个球面上的。这很正常，因为日常生活中最典型的网格就是地球上的经纬线。但不幸的是，地球的经纬线并不能组成上面定义的二维网格。因为按上面的定义，在网格上，往上下左右四个方向都能够无限地前进，而地球却有南北两极。
「吃豆子」游戏（ Pacman ）的地图符合上面的定义，因为从地图的左（右）边出去可以从右（左）边回来，从上（下）边出去可以从下（上）边回来。你可曾想过，这样的地图怎样才能实现？答案是环面：
photo-media
张文琦同学用欧拉示性数证明了满足上面定义的二维网格所「撑」起来的曲面，必定同胚于环面。通过计算机绘图观察，他认为二维网格经过迭代也收敛成环面。
三、用傅里叶变换分析迭代结果的形状
我用 Matlab 绘制得到的一次迭代结果如下图所示。我们看到，图形打「结」了，而且这个「结」经过了上千次迭代都没有消失。
photo-media
怎么回事呢？我决定再次拿出傅里叶变换这个大杀器。与多边形的情形一样，二维网格的迭代关于 x、y、z 三个坐标也是独立的。所以，我们可以只分析 x 坐标。
在分析多边形时，我们把各个点的 x 坐标看成了一个离散的、周期性重复的一维序列。在这里，我们可以把二维网格上各点的 x 坐标看成一个离散的、周期性重复的二维序列。而迭代操作，则是把这个二维序列与二维滤波器 [1/4, 1/4; 1/4, 1/4] 进行二维循环卷积。
这个滤波器的二维傅里叶变换的模如下：
photo-media
除去直流分量（左上角），模最大的两个值位于第 2 行第 1 列（以及左下角）和第 1 行第 2 列（以及右上角）。迭代过程中能留下来的也正是这两个分量，它们在原域中分别是 sin(2pii/n + P1)以及 sin(2pij/m + P2)（ P1、P2 为任意相位）。于是我们可以得到迭代足够多次后的网格表达式：
x(k)[i][j] ≈ Cx1(k) * sin(2pii/n + Px1(k)) + Cx2(k) * sin(2pij/m + Px2(k)) y(k)[i][j] ≈ Cy1(k) * sin(2pii/n + Py1(k)) + Cy2(k) * sin(2pij/m + Py2(k)) z(k)[i][j] ≈ Cz1(k) * sin(2pii/n + Pz1(k)) + Cz2(k) * sin(2pij/m + Pz2(k))
其中 Cx1、Cy1、Cz1、Cx2、Cy2、Cz2、Px1、Py1、Pz1、Px2、Py2、Pz2 由初始网格决定，在迭代过程中按滤波器的特性演变。
注意到上面三个式子中，如果仅保留右边的第一项，则可以得到一个三维空间中的椭圆，其中心位于原点；仅保留右边的第二项，也可以得到一个中心位于原点椭圆。于是极限网格就是由这两个椭圆「叠加」而成的。精确地说，是把一个椭圆 A 的中心放到另一个椭圆 B 上，然后沿着椭圆 B 平移椭圆 A 所生成的曲面。下面我画出了一次迭代过程，可以看到两个椭圆渐渐成形。曲面的成形则略慢一些。
photo-media
四、极限曲面的性质及与环面的区别
与「一个椭圆沿另一个椭圆平移得到的曲面」有关的研究似乎并不多。我在 Google 上搜索这种曲面的名称，只发现了 2013 年 6 月 arXiv 上的一篇文章，其中称这种曲面为 elliptic translational celestial。
首先我们来研究一下 elliptic translational celestial 这个名称。Elliptic 是「椭圆的」， translational 是「平移的」，这两个词比较好理解。celestial 我没有找到中文翻译，但这个词出现在「celestial body」（天体）中，也许来自于天文学。事实上，天文学中的确可以找到这种曲面，例如在地球的公转过程中，地球的赤道（椭圆 A ）沿地球公转轨道（椭圆 B ）平移形成的曲面就是 elliptic translational celestial （忽略地球自转引起的赤道自转），只是我不知道在天文学上这种曲面有没有研究价值。
Elliptic translational celestial （下文简称 ETC ）与环面看起来都是一个椭圆沿另一个椭圆运动所形成的曲面，那么二者有什么区别呢？区别在于，形成环面时椭圆 A 并不是在平移，而是在绕椭圆 B 的法线旋转。如果你想象一个钥匙链（椭圆 A ）套在一个呼啦圈（椭圆 B ）上，那么你沿着呼啦圈去撸钥匙链，钥匙链的运动大约就是旋转。那么在实际中，钥匙链能平移么？不行！如果真要平移，在呼啦圈上的两个位置，钥匙链必须穿过呼啦圈。正是由于这个原因，在 ETC 的两头会形成两个「结」。
photo-media
请看上图，这是椭圆 A （蓝色）沿椭圆 B （红色）平移所形成的 ETC。乍一看很像个环面，但请注意它的左右两端是打了结的，红色椭圆从这个结里露了出来。如果蓝色椭圆再「胖」一些，这个结就会更加明显：
photo-media
在这张图中，似乎「洞」（上下）与「结」（左右）的区别已经不那么明显了。对！如果蓝色的椭圆变得更大，大过了红色的椭圆，那么这个图就更容易被理解成红色椭圆沿着蓝色椭圆平移而形成。这时，上下的连通「洞」被截断成为两个「结」，而左右的两个「结」打通形成一个「洞」：
photo-media
总结一下 ETC 的性质，其实是这样的： 1) 两个椭圆各有一条法线，在这两条法线方向上会形成「洞」或「结」； 2) 存在第三个方向与两条法线都垂直。在这个方向上「直径」比较大的椭圆，它的法线方向上是一个「洞」；在这个方向上「直径」比较小的椭圆，它的法线方向上是两个
突然一下冒出来好多好多讨论这个东西的。其实这个东西相当的 trivial 啊，学过线性代数就能搞定了，搞不懂居然还有人写 paper 讨论这个东西...
首先，假设我们取了 n 个随机的点，把它们的坐标记做：
photo-media
然后，假如进行一次「取中点」操作以后，新的点的坐标记做：
photo-media
那么，显然有：
photo-media
可以把它们写成矩阵的形式：
photo-media
嗯。为了简单起见，把那个大方阵记做 M. 我们把取中点操作 k 次以后的坐标点记为：
photo-media
那么有：
photo-media
注意到横坐标和纵坐标其实是不相关的，也就是说：
photo-media
photo-media
那么要知道取了很多次中点以后，是什么情况，只需要知道 M 这个矩阵乘 k 次以后是什么样子的就行了。很自然，算 M 的特征值。
photo-media
这个矩阵的特征多项式恰好超级简单，就是
photo-media
所以写出所有的 n 个特征值：
photo-media
这些特征值在复平面上均匀地分布在圆心为 1/2 + 0i , 半径为 1/2 的圆上。特别注意到 1 总是特征值之一。
下面看一下特征向量：
photo-media
注意到：
photo-media
所以可以得到第 j 个特征向量的第 k 个分量和 k+1 分量的递推关系：
photo-media
最后，取所有特征向量的第一个分量为 1 ，就得到：
photo-media
好，既然有 n 个特征向量而无简并，则这些特征向量是完备基，也就是说任意一个 n 维列向量都可以写成这些特征向量线性组合，假设 r 是任意一个 n 维列向量，就有：
photo-media
那么由特征向量的定义，有：
photo-media
注意到，模最大的特征值就是 1 了，其他特征值的模都小于 1. 对于特征值为 1 的这个特征向量，它其实是：
photo-media
也就是说它的唯一作用就是使得整个图形在不断取中点的过程中，重心保持不变。
由于其他特征值的模都小于 1 ，所以整个图形随着你不断取中点肯定会慢慢缩成一个点(正好是它们的重心)，但是如果你每一取一次做一个归一化(normalization)，那它就能慢慢变成椭圆。原因在于模第二大的特征值。
模仅小于 1 的特征值有两个，从
photo-media
看出，模第二大的特征值是 j=2 的情形，以及其复共轭(因为 M 是个实矩阵，所以任一特征值的复共轭也是特征值)。为了下面讨论方便，我们给一些特殊的记号，把 j=2 的特征值记为：
photo-media
对应特征向量记为：
photo-media
注意到除了 1 以外，
对 photo-media
来说，其他的特征值得模都小于它们。
我们重新把任意向量 r 写成如下的形式：
photo-media
那么，有：
photo-media
这里打省略号，是因为其他特征值的模都比
photo-media
的小，所以在 k 足够大的时候与列出的这三项相比就趋于 0 了。
为了讨论方便，再记：
photo-media
photo-media
那么再代入
photo-media
就可以得到：
photo-media
大功告成，我
...
9300 用傅里叶变换玩转椭圆来源： 王赟 ^_^Maigo 的日志 这周人人网上有一张很火的照片，描述了这样一个有趣的数学现象：
在平面上任取 n 个点连成一个多边形（可以自交），然后取多边形各边的中点连成一个新多边形，不断重复上述操作，图形会收敛到一个椭圆。原图还附有一篇论文，用线性代数的方法证明了上述结论。
论文比较长，虽然排版很稀疏，但也长达 19 页。后来在人人网上看到这篇帖子，从信号处理的角度来解释椭圆是怎么形成的，让我豁然开朗。可惜这篇帖子对于傅里叶变换的部分并没有讲得很清楚，还引入了并无必要的高斯分布。所以我在这里用简洁的语言说明一下椭圆形成的道理，然后再附上一些好玩的拓展。
一、现象的严格表述
设初始时 n 个点的横坐标为 x(0)[0], x(0)[1], ..., x(0)[n-1] ，纵坐标为 y(0)[0], y(0)[1], ..., y(0)[n-1]。（注：为排版方便，我用小括号表示迭代次数，中括号表示顶点编号）
第 k 次操作如下： x(k)[i] <- ( x(k-1)[i] + x(k-1)[i+1] ) / 2, y(k)[i] <- ( y(k-1)[i] + y(k-1)[i+1] ) / 2 （对于 i=n-1 ，取 i+1=0 ）
按这种操作，其实点列不会收敛到椭圆，而是收缩到一个点，这个点是初始多边形各个点的平均值（或称重心，但注意这是顶点集的重心，不是多边形面积的重心）。为了避免图形收缩成一个点，需要对迭代过程做如下修正：
1) 初始时，通过平移将顶点集的重心移到原点； 2) 每次迭代后，在 x、y 方向上分别进行缩放，使得所有顶点的横坐标平方和为 1 ，纵坐标平方和也为 1。
这样修正后，图形依然不会收敛成一个椭圆（废话，它永远是个多边形），但迭代次数趋于无穷时，多边形的顶点会趋于排列在一个椭圆上。论文中还指出，这个椭圆的长轴与坐标轴成 45 度角。
二、用傅里叶变换解释椭圆的形成
首先注意到，迭代过程是对横、纵坐标独立地进行操作的，所以我们只要研究横坐标，对于纵坐标可以得到相同的结论。
我们把第 k 次迭代后得到的 n 个顶点的横坐标看作一个序列 x(k) = [ x(k)[0], x(k)[1], ..., x(k)[n-1] ] ，并想象对这个序列进行周期延拓，周期为 n。每次迭代操作相当于把延拓后的序列与序列 [1/2, 1/2] 进行卷积；在一个周期上，可以看成由 x(k)与序列 a = [1/2, 1/2, 0, ..., 0] （ n-2 个 0 ）进行循环卷积得到 x(k+1)。好了，原域的循环卷积在傅里叶变换域中就是乘法。我们用大写字母表示对应序列的傅里叶变换，则有
X(k+1) = A X(k) （这里的乘法表示对应元素相乘，下同）
于是就有
X(k) = A^k X(k)
所以，一切的奥妙就都蕴含在序列 a 的傅里叶变换 A 的性质里了。
学过信号处理的同学一定看出来了， a 是一个平滑滤波器，具有低通的性质。它的傅里叶变换为：
A[i] = 1/2 * (1 + exp(-i*2pi/n))
它的模是
abs(A[i]) = sqrt(1/2 + 1/2 * cos(i2pi/n)) = abs(cos(ipi/n))
画出来是这个样子（我取了 n=64 ）：
photo-media
嘿，还真是低通啊。直流分量没有变化，最高频分量（位于横轴正中央）一下子就被消灭干净了。
把这个低通滤波器反复施加于原始序列 x(0)，就可以想象到会发生什么了：直流分量始终保持不变，交流分量逐次衰减，越是高频衰减得越快。于是，在极限情况下， x 和 y 序列就都只剩下了直流分量。这意味着所有点的横、纵坐标都分别相同了，即图形收缩成一个点。这个直流分量的值正是初始多边形各点坐标的平均值，即重心。
我们现在加上修正 1。修正 1 的作用，就是去掉了序列中的直流分量。于是「山中无老虎，猴子称霸王」，最低频分量现在占了统治地位。尽管它也在衰减，但更高频分量衰减得更快。这时我们再看修正 2 ，它其实是让序列的总「能量」保持不变。在迭代次数足够多，高频分量衰减得可以忽略不计的时候，整个序列的能量就是最低频分量的能量，也就是说，在极限中被保留下来的就是恒定能量的最低频分量。最低频分量是一个周期为 n 的正弦函数，写成表达式就是：
x(k)[i] ≈ Cx(k) * sin(i2pi/n + Px(k)) y(k)[i] ≈ Cy(k) * sin(i2pi/n + Py(k))
（这里 Cx、Cy 表示振幅， Px、Py 表示相位）
修正 2 中说的让横、纵坐标的平方和均为 1 ，在迭代次数趋于无穷时，其实就是要求 x(k)[i] 和 y(k)[i] 的振幅相同。于是 x(k)[i] 和 y(k)[i] 是两个同频、同幅、不一定同相的正弦波。把它们看成关于 i*2pi/n 的参数方程，就不难发现它是一个椭圆，极限情况下多边形的各点就分布在这个椭圆上。由于两个正弦波振幅相同，这个椭圆的长轴与坐标轴成 45 度角。这个椭圆的形状由相位差 Px(k) - Py(k)决定。由于每次迭代对 x 和 y 相位的影响是相同的，所以这个相位差就是初始序列 x(0)[i]、y(0)[i] 中最低频分量的相位差。
* 当这个相位差为 0 度时，椭圆退化成 y=x 上的一条线段； * 当这个相位差在 0~90 度之间时，椭圆长轴沿 y=x 方向，相位差越大，椭圆越「胖」； * 当这个相位差恰为 90 度时，椭圆变成正圆； * 当这个相位差在 90~180 度之间时，椭圆长轴沿 y=-x 方向，相位差越大，椭圆越「瘦」； * 当这个相位差为 180 度时，椭圆退化成 y=-x 上的一条线段。
三、一些好玩的拓展
- 高维情况
由于各维的变换完全独立，在高维中，图形也会收敛成一个椭圆，其长轴的方向向量为(±1, ±1, ..., ±1)，说白了就是跟各个坐标轴的夹角都相等。有人已经做出了三维情况下的动画。
- 如果每次迭代不是取线段中点，而是取相邻三个顶点的重心会怎样？
此时 a = [1/3, 1/3, 1/3] （后面的 0 省略），其傅里叶变换的模如下：
photo-media
所以，「猴子」还是最低频分量，结论不变，图形仍收敛于椭圆。
- 如果每次迭代是取相隔一个顶点的两个顶点连线的中点为新的顶点，会怎样？
此时 a = [1/2, 0, 1/2] ），其傅里叶变换的模如下：
photo-media
好玩的事情来了！现在这是一个带阻滤波器。当 n 是偶数的时候，有一个最高频分量不会衰减，于是原来的最低频分量就怂了。这个最高频分量的周期是 2 个点啊！于是，原多边形中所有编号为奇数的点会收缩成一个点，编号
Q. E. D.
