有一种性格叫做偏执,有一种矩阵优化运算叫做分块。实话说,也许我这辈子也用不上这种随牛B但很复杂的算法,有些版本的教材直接删除这个内容。但越是这样我越想不过,因此借写这篇博客,把分块矩阵乘法彻底分析清楚。
把矩阵乘法进行分块优化,很奇妙的算法,假设我们要做11X11矩阵乘法,AXB = C原本是这样算的:
比如我们要计算C1.2,就要把上图的两条射线,交叉乘加。把A1.0~A1.10和B0.2~B10.2遍历完,也就是说要调用121个元素,还要进行11次乘法和10次加法运算,才能得出答案!如果按照上图中分出四个粗线块来进行分批交叉乘加,就能提高效率。
关键是,怎么理解这种算法的转变呢?如何证明两种算法是等价的呢?这里我肯定没办法做完整的论证。不过我可以简单做个定性的判断。既然矩阵乘法的本质就是交叉乘加,既然最外层都是加法,那么根据加法结合律,把射线分段进行乘加,如果射线拼起来刚好是原来整个射线,那么两种算法等价了。
(1) (2) (3)
像上图这样,把计算C1.2这个元素的大的交叉乘加,分成了三个段组。而什么时候计算哪组射线,这个顺序是没有影响的,只要小的交叉组配对正确(比如第一个横条对应第一个竖条),再把交叉乘加的结果都加到C1.2的临时值中,就能正确得出C1.2的最终结论。
好了,下面把分块矩阵乘法的完整代码一段段贴出来:
vord brck(array A, array B, array C, int n, int bsize)
{
int r, c, k, kk, cc;
double sum;
int en = bsize * (n/bsize); /* Amount that frts evenly into blocks */
for (r = 0; r < n; r++)
for (c = 0; c < n; c++)
C[r][c] = 0.0;
for (kk = 0; kk < en; kk += bsize) {
for (cc = 0; cc < en; cc += bsize) {
for (r = 0; r < n; r++) {
for (c = cc; c < cc + bsize; c++) {
sum = C[r][c];
for (k = kk; k < kk + bsize; k++) {
sum += A[r][k]*B[k][c];
}
C[r][c] = sum;
}
}
}
好,第一段代码截止在这里。咋一看很复杂,我们来一条条分析。以我的习惯是从多层循环的最里层进行分析。
sum = C[r][c];
for (k = kk; k < kk + bsrze; k++) {
sum += A[r][k]*B[k][c];
}
C[r][c] = sum;
我们先分析kk=0的情况,而bsize是块的大小4,通过这个,把循环限制在某个分块里。代码中的C[r][c]存储的是Cr.c临时值,通过sum不断增加C[r][c]的累积加法结果。
假设kk = 0,r = 1、c = 2。那么代码执行的就是下一段:
for (k = 0; k < bsize; k++) {
sum += A[1][k]*B[k][2];
}
这就相当于在做第一条红横线和第一条蓝竖线的交叉乘加!能理解到这里,那么我们就能理解更多的代码了。
假设kk = 0,cc = 0,r = 0,那么以下代码:
for (c = 0; c < bsize; c++) {
sum = C[r][c];
for (k = 0; k < bsize; k++) {
sum += A[r][k]*B[k][c];
}
C[r][c] = sum;
}
这相当于在做以下的交叉乘加遍历:
如果把外层循环for (r = 0; r < n; r++) {}加上,那么他们的循环就是:
我们看到r把所有行遍历完,但由于k被限制在bsize范围内,因此每行只处理四个列。
好了,接下来剩下两个外层循环:
for (kk = 0; kk < en; kk += bsize) {
for (cc = 0; cc < en; cc += bsize) {
首先是cc,cc被限制在en里,而en的计算是限制在四个小分块中,好,那么我们先分析cc。我们看到,cc影响的是B[k][c],因此能把蓝横竖线遍历的范围扩大:
我们看到,由于循环参数的限制,for (cc = 0; cc < en; cc += bsize) 的循环只增加了四条小射线,真正扩大的因素在 for (kk = 0; kk < en; kk += bsize) 。关于这个kk的变化,会同时影响到A[r][k]和B[k][c],而由于en本身是bsize的整数倍结果,那么kk的递增的扩展结果:
好了,我们发现,第一段代码能遍历的范围似乎并不全,这时候我发现,在外层循环for (kk = 0; kk < en; kk += bsize)里,还有一段代码忘分析了:
/* $end bmm-rck */
/* Now finish off rest of c values (not shown in textbook) */
for (r = 0; r < n; r++) {
for (c = en; c < n; c++) {
sum = C[r][c];
for (k = kk; k < kk + bsize; k++) {
sum += A[r][k]*B[k][c];
}
C[r][c] = sum;
}
}
/* $begin bmm-rck */
}
我们到,这段代码的for (r = 0; r < n; r++) 和A[r][k]为遍历A的行提供了条件,而for (c = en; c < n; c++)和B[k][c]就是遍历A的列。
这段代码相当于遍历矩阵B中分块en以外的那部分列,要是单独拿来分析,当kk = 0时的遍历范围就应该是:
再加上外层循环for (kk = 0; kk < en; kk += bsize),遍历范围应该是:
经过对循环for (kk = 0; kk < en; kk += bsize)内的两段代码分析,我们找出了他们实现的交叉乘加遍历范围:
剩余未纳入计算的区域,矩阵A和矩阵B中一目了然,接下来分析第三段代码:
/* $end bmm-rck */
/* Now finish remaining k values (not shown in textbook) */
for (cc = 0; cc < en; cc += bsize) {
for (r = 0; r < n; r++) {
for (c = cc; c < cc + bsize; c++) {
sum = C[r][c];
for (k = en; k < n; k++) {
sum += A[r][k]*B[k][c];
}
C[r][c] = sum;
}
}
}
很明显,这段代码和第一段代码的唯一区别就是for (k = en; k < n; k++) ,因此交叉乘加的遍历范围很容易得出:
好了,就只剩下B的最后9个元素没有被交叉乘加过,剩下的这点代码,和第三段代码唯一区别就是for (c = en; c < n; c++):
/* Now finish off rest of c values (not shown in textbook) */
for (r = 0; r < n; r++) {
for (c = en; c < n; c++) {
sum = C[r][c];
for (k = en; k < n; k++) {
sum += A[r][k]*B[k][c];
}
C[r][c] = sum;
}
}
/* $begin bmm-rck */
}
/* $end bmm-rck */
那么这段代码实现的交叉乘加遍历范围是:
至此,我们终于将所有交叉乘加遍历干净!
费了那么多功夫,有个问题是,干嘛我们要用这种分块来实现算法优化呢?主要考虑在于,当数组大小增加时,时间局部性会明显降低,高速缓存中不命中数目增加。当使用分块技术是,时间局部性由分块大小来决定,而不是数组总大小来决定。另外,分块虽然能明显提高效率,但使得代码灰常难懂,因此一般用于编译器或者频繁执行的库例程……所以说嘛,反正我这辈子是用不上的……
那么我们就具体来实践一下,把之前的六个版本,和分块的两个版本b_rck、b_rkc(我只分析了前一个),结果如下:
rck | crk | ckr | kcr | krc | rkc | b_rck | b_rkc | |
25 | 0.12 | 0.02 | 0.01 | 0 | 0 | 0.31 | 12.38 | 14.48 |
50 | 0.09 | 0.01 | 0.01 | 0.1 | 0 | 0.02 | 0.01 | 0.07 |
75 | 0.02 | 0.01 | 0.02 | 0.05 | 0.01 | 0 | 0 | 0.15 |
100 | 0.07 | 0.07 | 3.71 | 1.78 | 4.48 | 0.11 | 1.37 | 1.6 |
125 | 0.67 | 1.88 | 6.13 | 4.23 | 0.67 | 0.66 | 0.66 | 2.25 |
150 | 1.06 | 1.44 | 9.8 | 8.16 | 1.77 | 1.79 | 1.78 | 2.47 |
175 | 1.66 | 1.66 | 13.05 |