优化程序中的数据：从数组到代数

前言

我们往往都希望优化我们的程序，使之达到一个更好的效果，程序优化的一个重点就是速度，加快速度的一个好办法就是使用并行技术，但是，并行时我们要考虑必须串行执行的任务，也就是有依赖关系的任务，任务中的重点往往是具体的数据，这些任务中的数据通常具有局部性和关联性。

而数据中数组具有代表性，现在，让笔者从数组开始，谈谈程序数据的优化。

从数据的存储内存开始

我们都知道计算机的基本内存结构如下：

而内存的结构又可以继续划分：

虚拟内存是一个很伟大的发明，它借助内存管理单元（MMU），并利用分页机制将磁盘的一部分模拟为内存使用。它允许计算机使用硬盘空间来扩展实际的物理内存。这使得操作系统能够运行超过实际物理内存容量的程序。

而我们重点关注的地方在cache这里。

cache命中率

cache会从更低一级的内存结构中搬数据，如果数据访问是局部性很强（如访问同一数据块多次），则缓存命中率会较高，如果不命中，那么计算机会跑到下一级内存中寻找数据，这样程序运行效率就会非常低。

优化

得知了这一点后，我们可以考虑改善我们的程序写法了，以数组操作为例：

for(int i = 0; i<= 2; i++){for(int j = i; j<= 2; j++)Z[j][i] = 0;
}

在C语言中，二维数组的内存分布通常是按行优先（Row-major order）存储的，这意味着数组的行是连续存储在内存中的。具体来说，对于一个二维数组 Z，其内存布局是按以下方式排列的：

二维数组的内存分布

假设我们有一个二维数组 Z，其大小为 m 行 n 列。数组元素在内存中的排列顺序如下：

Z[0][0], Z[0][1], ..., Z[0][n-1], Z[1][0], Z[1][1], ..., Z[1][n-1], ..., Z[m-1][0], ..., Z[m-1][n-1]

每一行的元素是连续存储的，然后依次存储下一行的元素。

那么，优化后的遍历方法如下：

for (int j = 0; j <= 2; j++) {for (int i = 0; i <= j; i++) {Z[j][i] = 0;}
}

上面的优化方法相信大家都能琢磨出来，但是，如果稍微改一下呢？

for(int i = 0; i<= 5; i++){for(int j = i; j<= 7; j++)Z[j][i] = 0;
}

按照原程序的遍历：

┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
│ 0,0 │ 1,0 │ 2,0 │ 3,0 │ 4,0 │ 5,0 │ 6,0 │ 7,0 │
│     │ 1,1 │ 2,1 │ 3,1 │ 4,1 │ 5,1 │ 6,1 │ 7,1 │
│     │     │ 2,2 │ 3,2 │ 4,2 │ 5,2 │ 6,2 │ 7,2 │
│     │     │     │ 3,3 │ 4,3 │ 5,3 │ 6,3 │ 7,3 │
│     │     │     │     │ 4,4 │ 5,4 │ 6,4 │ 7,4 │
│     │     │     │     │     │ 5,5 │ 6,5 │ 7,5 │
│     │     │     │     │     │     │ 6,6 │ 7,6 │
│     │     │     │     │     │     │     │ 7,7 │
└─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘

更好的遍历方法：

┌─────┬─────┬─────┬─────┬─────┬─────
│ 0,0 │ 										
│ 1,0 │ 1,1 │  									
│ 2,0 │ 2,1 │ 2,2 │  
│ 3,0 │ 3,1 │ 3,2 │ 3,3 │ 
│ 4,0 │ 4,1 │ 4,2 │ 4,3 │ 4,4 │ 
│ 5,0 │ 5,1 │ 5,2 │ 5,3 │ 5,4 │ 5,5 
| 6,0 | 6,1 | 6,2 | 6,3 | 6,4 | 6,5 
| 7,0 | 7,1 | 7,2 | 7,3 | 7,4 | 7,5 
└─────┴─────┴─────┴─────┴─────┴─────

局部性更好的程序如下，此时想要一眼看出来这样写就有点困难了，那我们要怎么推导数组的遍历式呢：

for (int j = 0; j <= 7; j++) {for (int i = 0; i <= (j < 5 ? j : 5); i++) {Z[j][i] = 0;}
}

引入线性代数

我们先看看各种值的范围：

i的范围: i>=0, i<=5
j的范围: j>=i, j<=7

尝试把它们写成线性方程：

1*i + 0*j + 0 >= 0
-1*i + 0*j + 5 >= 0-1*i + j + 0 >= 0
0*i + -1*j + 7 >= 0

矩阵如下：

|  1  0 |   | i |   >=   | 0 |
| -1  0 | * | j |   >=   | -5 |
| -1  1 |           >=   | 0 |
|  0 -1 |           >=   | -7 |

现在我们得到了矩阵，我们可以进一步得到多面体，先回顾一下矩阵与多面体的关系：

线性约束表示多面体

多面体可以通过一组线性不等式来定义，这些不等式可以表示为矩阵和向量的形式。例如，对于一个包含 n个变量的多面体，可以用一个 m×n 的矩阵A和一个m维的向量 b来表示：

Ax <= b

其中，x是变量向量，约束条件定义了多面体的边界。

顶点表示

多面体的顶点可以通过求解线性方程组（通常涉及矩阵的逆或者伪逆）来获得。这些顶点是满足约束条件的解。

矩阵操作多面体

线性变换

通过矩阵乘法，可以对多面体进行线性变换（如旋转、缩放、平移等）。例如，如果矩阵M描述了一个线性变换，那么多面体中的每一个点 x在变换后的新位置可以表示为Mx。

仿射变换

仿射变换是线性变换的推广，包括线性变换和平移。可以用如下形式表示：

y=Mx+t

其中，MM 是线性变换矩阵，t是平移向量。

好吧，其实矩阵和多面体与接下来要讲的算法也没多大关系，笔者只是想说明如何从不等式推导到线性代数并扩展到多面体和高维空间体的。

使用Fourier-Motzkin算法

Fourier-Motzkin算法是一种经常在多面体中用于求解线性不等式系统的消去算法，概括如下：

选择消去变量： 选择一个变量 xi作为消去变量。

分类不等式： 将所有不等式分为三类：

• 包含 Xi的不等式，且 Xi的系数为正。
• 包含 Xi的不等式，且 Xi的系数为负。
• 不包含 Xi的不等式。

生成新不等式： 通过将第一类不等式和第二类不等式配对，消去 Xi

组合不等式： 将生成的新不等式与不包含 xi的不等式组合，得到一个新的线性不等式系统。

重复步骤： 对新的线性不等式系统重复上述步骤，直到所有变量都被消去。

应用该算法，我们重新得到范围：

0<=j, 0 <=5
j<=7

那么i和j的范围如下：

L(i):0
U(i):5,jL(i):0
U(J):7

有了这个范围，我们可以得到：

for (int j = 0; j <= 7; j++) {for (int i = 0; i <= min(5,j); i++) {Z[j][i] = 0;}
}

也就是：

for (int j = 0; j <= 7; j++) {for (int i = 0; i <= (j < 5 ? j : 5); i++) {Z[j][i] = 0;}
}

总结

从程序中的优化出发，由程序存储引出cache，再由cache命中率引出数据局部性的重要性，为了提高数据局部性，必须改变循环遍历方法。为了改变循环遍历方法，由不等式引出线性代数，再由线性代数引出多面体，最后使用算法计算约束，得到具有良好局部性的程序。
其实没啥好总结的，只写了一小段，还没写完开头呢，不过先更到这，该上床睡觉了。