循环矢量化

本页目录

1．数组运算矢量化方法

MATLAB是一个矩阵语言，是为矢量和矩阵操作设计的，可以通过矢量化方法来加速M文件的运行。矢量化指的是将for循环和while循环转换为等价的矢量或矩阵操作。例如：下面以0.01为间隔计算0到10之间1001个数的正弦值。

code.matlab 
i = 0;
for t = 0:.01:10
    i = i+1;
    y(i) = sin(t);
end

矢量化以后，为：

code.matlab 
t = 0:.01:10;
y = sin(t);

矢量化以后比未矢量化时要快很多，可以用tic和toc函数进行测试。

数组运算对数组中的每个元素进行相同的运算。计算1000个三角形的面积，用循环结构进行计算：

code.matlab 
for i=1:1000
  S(i)=B(i)*H(i)/2
end

其中，B(i)和H(i)分别为三角形的底和高的长度。

矢量化后的计算为：

code.matlab 
S=B.*H/2

2．逻辑数组运算

可以对数组进行逻辑运算。比如下面的数组中，元素0表示女性，元素1表示男性。

code.matlab 
>> gender=[1 0 0 1 1 1 0 1 0 1];
>> gender==0
ans =
     0     1     1     0     0     0     1     0     1     0

下面的math数组表示这些学生的数学成绩，要求显示其中女生的成绩。

code.matlab 
>> math=[89 90 78 91 69 83 80 76 89 95];
>> mathf=math(gender==0)
mathf =
    90    78    80    89

下面显示成绩大于85分的男生的成绩。

code.matlab 
>> mathf85=math((math>85) & (gender==1))
mathf85 =
89  91    95

3．矩阵运算

下面根据给定的x和y值，计算二维函数曲面网格节点上的z值。假设函数为z=sin(x)*cos(y)。

采用循环结构进行计算。

code.matlab 
i=0;
j=0;
for m=-3:0.1:3
  i=i+1;
  j=0;
  for n=-3:0.1:3
    j=j+1;
    X(i,j)=m;
    Y(i,j)=n;
    Z(i,j)=sin(X(i,j))*cos(Y(i,j));
  end
end

矢量化以后，表示为：

code.matlab 
x=-3:0.1:3;
y=-3:0.1:3;
[X,Y]=meshgrid(x,y);
Z=sin(X).*cos(Y);

meshgrid函数是一个实现了矢量化的函数，它简化了上面两层嵌套的for循环。可以用下面的语句编辑该函数的代码。

code.matlab 
edit meshgrid

下面用矢量化的方法构造矩阵。

code.matlab 
>> v=1:4;
>> M=repmat(v,4,1)
M =
     1     2     3     4
     1     2     3     4
     1     2     3     4
     1     2     3     4

repmat函数是一个实现了矢量化的函数。

4．排序和计数运算

可以用矢量化的方法实现排序和计数等操作。

code.matlab 
>> x=[2 4 3 6 4 1 9 4 8 2];
>> x=sort(x);
x =
     1     2     2     3     4     4     4     6     8     9

计算x中值为4的元素的个数。

code.matlab 
>> sum4=sum(x==4)
sum4 =
     3

5．使用了矢量化的函数

MATLAB中，有些函数内部已经采用了矢量化处理，因而运行效率比较高。这些函数的函数名如表1-1所示。

表1-1 MATLAB中使用了矢量化的函数

all diff ipermute permute reshape squeeze
any find logical prod shiftdim sub2ind
cumsum ind2sub ndgrid repmat sort sum
meshgrid

下面以repmat函数为例，具体介绍它内部是如何进行矢量化处理的。该函数有3个输入变量：数组A、行维数M和列维数N，通过将数组A做M×N的“叠置”操作来生成新的矩阵，如下所示：

code.matlab 
>> A = [1 2 3; 4 5 6];
>> B = repmat(A,2,3)
B =
    1    2    3    1    2    3    1    2    3
    4    5    6    4    5    6    4    5    6
    1    2    3    1    2    3    1    2    3
    4    5    6    4    5    6    4    5    6

repmat函数通过矢量化创建编码，这些编码把元素放在输出数组中。该函数的代码如下所示：

code.matlab 
function B = repmat(A,M,N)
if nargin < 2
   error('Requires at least 2 inputs.')
elseif nargin == 2
   N = M;
end
% 第1步获取行和列的大小
[m,n] = size(A);
% 第2步生成编号从1到行或列的大小的矢量
mind = (1:m)';
nind = (1:n)';
% 第3步根据上面的矢量创建编号矩阵
mind = mind(:,ones(1,M));
nind = nind(:,ones(1,N));
% 第4步创建输出数组
B = A(mind,nind);

可见，repmat函数通过4个步骤来生成新矩阵。

第1步获取输入数组行和列的大小。

第2步创建两个列矢量，mind包含从1到A的行数的整数，nind包含从1到A的列数的整数。

第3步通过一个MATLAB矢量化技巧把单列数据复制到任意多列。代码为：

B = A(:,ones(1,n_cols))

其中，n_cols是最终生成的矩阵的列数。

第4步通过数组索引创建输出数组。行号数组mind的每个元素与列号数组mind的每个元素相匹配，通过下面的过程来实现：

（1）mind的第1个元素，即行号，与nind的每个元素相匹配。MATLAB按列序在nind矩阵中移动，所以mind(1,1)首先从nind(1,1)开始，然后是nind(2,1)。其结果是输出数组的第1行被填满。

（2）对mind中的每个元素重复上面的过程，得到输出数组。