众所周知，向量化编程在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' http://projecteuler.net/images/symbol_maps.gif精华 · 高效向量化函数" /> 11
11 is 'two ones' http://projecteuler.net/images/symbol_maps.gif精华 · 高效向量化函数" /> 21
21 is 'one two and one one' http://projecteuler.net/images/symbol_maps.gif精华 · 高效向量化函数" /> 1211
1211 is 'one one, one two and two ones' http://projecteuler.net/images/symbol_maps.gif精华 · 高效向量化函数" /> 111221
111221 is 'three ones, two twos and one
one' http://projecteuler.net/images/symbol_maps.gif精华 · 高效向量化函数" /> 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 =
1012.
Give your answer modulo 230 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函数尽管运算时间较长，却可以使代码简化。在一些特殊情况下，它们可以达到非常高效的运算效果。
加载中…
(2014-04-13 18:01:27) 