【讨论】大型分块矩阵求和算法效率的比较

winner245

问题：现有一个M*M的大矩阵，欲对其每个分块子矩阵N*N求和，其中，b = M/N 是一个整数，所得和矩阵为 b*b。下面给出几种实现算法，并对其效率做比较，希望能抛砖引玉，找到更高效的算法。

现象：分别用 accumarray 法和普通的 "sum + reshape" 法计算了这个分块矩阵求和，发现普通的算法 "sum + reshape" 效率比 accumarray 高 20~30 倍左右，这与我们通常认为的 accumarray  效率很高似乎不符。

效率比较：对于 accumarray 函数，需要事先生成 subs，由于 subs 可能为线性索引或下标索引，accumarray 法又细分为两类方法：线性索引+accumarray ，下标索引+accumarray 。为了充分说明 accumarray 法效率更低，在测试代码效率时，没有计入生成 subs 时的额外开销（故若考虑subs生成的开销，accumarray 的效率还会更低）。另外，分别使用了 tic/toc 法和 timeit 法测试比较，所得结果基本吻合。

平台：MATLAB 2014a， i7-3630QM，64bit Win 8

基于 tic/toc 的效率比较（程序在运行多遍下的测试结果）：

1. M = 5000; N = 100; b = M/N;
   
2. a = randi(10,M);
   
3. subs1 = reshape(kron(reshape(1:b^2,b,b),ones(N)),[],1);
   
4. subs2 = reshape(repmat(kron(reshape(fullfact([b b]),b,[]),ones(N,1)),N,1),[],2);
   
5. tic, sum1 = reshape(accumarray(subs1,a(:)),b,b); t1 = toc    % 线性索引+accumarray
   
6. tic, sum2 = accumarray(subs2,a(:)); t2 = toc    % 下标索引+accumarray
   
7. tic, sum3 =reshape(sum(sum(reshape(a,N,b,N,b),1),3),b,b); t3 = toc  % 普通算法 "sum + reshape"
   
8. [t1/t3, t2/t3]
   
9. [isequal(sum1,sum3), isequal(sum2,sum3)]

复制代码

t1 =
    0.2542

t2 =
    0.3327

t3 =
    0.0126

ans =
   20.1760   26.4092

ans =
     1     1

基于 timeit 的效率比较：

1. M = 5000; N = 100; b = M/N;
   
2. a = randi(10,M);
   
3. subs1 = reshape(kron(reshape(1:b^2,b,b),ones(N)),[],1);
   
4. subs2 = reshape(repmat(kron(reshape(fullfact([b b]),b,[]),ones(N,1)),N,1),[],2);
   
5. sum1 = @() reshape(accumarray(subs1,a(:)),b,b);  % 线性索引+accumarray
   
6. sum2 = @() accumarray(subs2,a(:));   % 下标索引+accumarray
   
7. sum3 =@() reshape(sum(sum(reshape(a,N,b,N,b),1),3),b,b);   % 普通算法 "sum + reshape"
   
8. t1 = timeit(sum1)
   
9. t2 = timeit(sum2)
   
10. t3 = timeit(sum3)
    
11. [t1/t3, t2/t3]
    
12. [isequal(sum2(),sum3()) isequal(sum2(),sum3())]

复制代码

t1 =
    0.2516

t2 =
    0.3345

t3 =
    0.0095

ans =
   26.5604   35.3143

ans =
     1     1

现象/讨论：

1. 可以看到即使在不计入 subs 生成开销的情况下，对于大型分块矩阵的求和，accumarray 的效率依然不及普通算法 "sum + reshape"。且效率相差约30倍左右。其实在求这个问题之初，之所以优先考虑 accumarray 函数，是因为 accumarray 的名字本身就是求和的意思，而且，通常认为 accumarray  效率很高，此处的效率比较结果还是有些出乎意料。不过，在得到效率比较结果后反观普通算法，只涉及了几次 reshape 和 sum，没有任何复杂函数，也许正因为如此，普通算法才更高效。

2. 两种 accumarray 方法里，基于线性索引的算法效率略高，说明在使用 accumarray 时，应优先考虑使用线性索引 （当然，如果生成线性索引的额外开销大于下标索引，则效率的比较还要考虑subs生成的开销）。

rocwoods

刘兄这个问题不错。手机上回复不方便，简单说下。accumarray的高效个人认为体现在一些计算复杂度高的场合，因为它本身的算法可能是有针对性的。而分块求和这种直觉上很容易想到算法的低复杂度的问题如果硬要用accumarray来解，就有点用牛刀杀鸡，大炮打蚊子的感觉了。刘兄给的第三种方法应该是最高效的之一了，如果对效率要求不是十分苛刻，也可以用blockproc，简单明了。

winner245

rocwoods 发表于 2014-9-15 18:43
刘兄这个问题不错。手机上回复不方便，简单说下。accumarray的高效个人认为体现在一些计算复杂度高的场合， ...

吴兄说得太好了，用 accumarray 来解低复杂度问题犹如牛刀杀鸡、大炮打蚊子了，一语惊醒梦中人啊！ blockproc 确实简明，一句话搞定 blockproc(a,[N N],@(block_struct) sum(block_struct.data(:)));