MATLAB 精华 · 高效向量化函数

2014-04-13 18:01阅读：

http://blog.sina.cn/dpool/blog/u/1560732217

众所周知，向量化编程在MATLAB中可以极大的提高程序运算效率。这其中的基本技巧包括： 1. 优化内存：预分配矩阵，避免动态变化矩阵大小等。
2. 向量化处理循环：使用尽可能多的向量化函数，尽量使用列向量。
3. 高效使用函数：内建函数优先级最高，其次是函数文件的p文件，最后才是函数句柄及内联函数。
即便for循环很多时候被要求尽量避免，MATLAB对其的优化已经使得for循环在很多小运算量的情况下仍然是效率最高的做法。除此之外，在矩阵运算中，高效的运用向量化函数可以极大的提高计算效率。以Project Euler的题目419为例，本文测试arrayfun, bsxfun及for循环的效率比较。
Project Euler 419：
The look and say sequence goes 1, 11, 21, 1211, 111221, 312211, 13112221, 1113213211, ...
The sequence starts with 1 and all other members are obtained by describing the previous member in terms of consecutive digits.
It helps to do this out loud:
1 is 'one one'

精华 · 高效向量化函数' /> 11
11 is 'two ones'

21
21 is 'one two and one one'

1211
1211 is 'one one, one two and two ones'

111221
111221 is 'three ones, two twos and one one'

312211
...
Define A(n), B(n) and C(n) as the number of ones, twos and threes in the n'th element of the sequence respectively.
One can verify that A(40) = 31254, B(40) = 20259 and C(40) = 11625.
Find A(n), B(n) and C(n) for n = 10¹².
Give your answer modulo 2³⁰ and separate your values for A, B and C by a comma.
E.g. for n = 40 the answer would be 31254,20259,11625
这个外观序列最初由Conway提出，并具备92个子序列。这些子序列可以进化成其他的子序列，因此整个序列的生成过程可以通过这个92*92大小的变化矩阵的乘法来实现。这其中涉及到的矩阵乘法运算次数仅仅由n来决定。对于题目的n，只要计算40次矩阵乘方即可。由于只需要计算结果的模，乘方过程中涉及到的最大二进制位数是60位。这样涉及到的一个主要问题就是，在MATLAB中，常用的计算数据类型是double型，而double的精度只可以达到2的53次方，因此只能采用uint64数据类型。但是MATLAB却仅仅支持uint64的基本运算，对于矩阵乘法，同余等运算均不支持。

底层的高效矩阵乘法算法的实现是基于硬件架构的，其指令包括各种hardware based优化。MATLAB对这些基本指令函数库进行了优化打包，形成方便的调用接口。因此，对于MATLAB不支持运算的数据类型，只有通过直接的计算方法。两个矩阵A*B=C的实现是分别对A和B的行列相乘得到C矩阵的每个元素。在MATLAB中，这可以通过对A的n行进行循环进行：
for line = 1:n
Atemp = repmat(A(line, :)', [1 n]); % Repeat each line to matrix
Ctemp = Atemp .* B;
C(line, :) = sum(Ctemp, 'native'); % Sum operation for native data type
end
这里使用点乘是因为MATLAB对所有的数据类型均支持点乘运算。

进一步的改进是使用bsxfun来替代repmat函数，因为repmat函数会带来内存复制的消耗。而bsxfun中的复制矩阵是采用高度优化的循环。bsxfun(@fun, a, b)的作用是对矩阵a和b进行fun运算，这个运算是对两个矩阵相应的元素进行。a和b可以都是或者一个是向量。MATLAB会自动复制向量成矩阵，和repmat效果一样，来进行矩阵的fun运算。因此代码可以化简为：
for line = 1:n
Ctemp = bsxfun(@times, A(line, :)', B); % Deal every line in A
C(line, :) = sum(Ctemp, 'native');
end
另外一种类似的函数arrayfun可以对矩阵批量的进行运算。通常arrayfun函数的使用是为了代码的简洁化，因为它需要调用大量的函数句柄。arrayfun(@fun, S)可以对S中的所有元素进行fun运算：
Acell = num2cell(A, 2); % Generate line cells of A
fun = @(x) sum(bsxfun(@times, x{:}', B), 'native');
C = cell2mat(arrayfun(fun, Acell, 'UniformOutput', false))
首先对A矩阵的所有行处理成cell结构，然后arrayfun可以对cell结构中的所有元素（A的每一行）进行处理，最后得到结果矩阵C。
此外，还有一种图像处理中常用的块处理函数blockproc也可以实现上述运算：
fun = @(blk) sum(bsxfun(@times, blk.data', B), 'native');
C = blockproc(A, [1 n], fun); % Process every line of A
这个函数的内部处理机制实际上是采用循环，对图像块进行处理。在图像处理中这种方法非常方便并且高效，因为很多场合下图像块间的处理需要有重合部分。
对于题目419中的92*92的矩阵进行39次乘法运算，在其中加入bitand函数进行uint64数据类型的同余运算，四种方法的运算时间为：
方法调用次数时间(s)
repmat 3588 0.34
bsxfun 3588 0.21
arrayfun 39 0.30
blockproc 39 0.48

可以看出，MATLAB对于基础for循环运算进行了非常好的优化，在对于小计算量的情况下，for循环仍然是最好的选择。然而，在这其中，尽可能多的向量化运算可以极大的提高运算效率，如使用bsxfun代替repmat。arrayfun和blockproc函数尽管运算时间较长，却可以使代码简化。在一些特殊情况下，它们可以达到非常高效的运算效果。

举报/Report

我的更多文章

下载客户端阅读体验更佳