向量化
向量化的应用
MATLAB® 针对涉及矩阵和向量的运算进行了优化。修正基于循环且面向标量的代码以使用 MATLAB 矩阵和向量运算的过程称为向量化。由于多种原因,代码的向量化是可取的:
外观:向量化的数学代码在外观上更像教科书中的数学表达式,从而使代码更便于理解。
不容易出错:由于不带循环,向量化的代码通常更短。代码行的减少意味着引入编程错误的机会也减少了。
性能:向量化代码的运行速度通常比相应的、包含循环的代码更快。
实现代码向量化以进行常规计算
以下代码计算 1,001 个从 0 到 10 之内的值的正弦:
i = 0; for t = 0:.01:10 i = i + 1; y(i) = sin(t); end
这是同一代码的向量化版本:
t = 0:.01:10; y = sin(t);
第二个代码示例的执行速度通常比第一个快,是 MATLAB 的更高效的用法。通过创建包含所示代码的脚本测试您系统上的执行速度,然后使用 tic
和 toc
函数测量其执行时间。
针对特定任务实行代码向量化
此代码计算某向量每五个元素的累加和:
x = 1:10000; ylength = (length(x) - mod(length(x),5))/5; y(1:ylength) = 0; for n= 5:5:length(x) y(n/5) = sum(x(1:n)); end
使用向量化,您可以编写一个更精确的 MATLAB 过程。以下代码显示了完成此项任务的一种方法:
x = 1:10000; xsums = cumsum(x); y = xsums(5:5:length(x));
数组运算
数组运算符对数据集中的所有元素执行相同的运算。这些类型的运算用于重复计算。例如,假设您通过记录各圆锥体的直径 (D
) 和高度 (H
) 来收集这些圆锥体的体积 (V
)。如果您只收集一个圆锥体的信息,则可以只计算该圆锥体的体积:
V = 1/12*pi*(D^2)*H;
现在,收集 10,000 个圆锥体的信息。向量 D
和 H
均包含 10,000 个元素,并且您需要计算 10,000 个体积。在绝大多数编程语言中,您需要设置类似于以下 MATLAB 代码的循环:
for n = 1:10000 V(n) = 1/12*pi*(D(n)^2)*H(n); end
借助 MATLAB,您可以使用与标量情形类似的语法对每个向量元素执行计算:
% Vectorized Calculation
V = 1/12*pi*(D.^2).*H;
注意
在运算符 *
、/
和 ^
之前放置句点 (.
),将其转换为数组运算符。
借助数组运算符,您还能够组合不同维度的矩阵。这样自动扩展大小为 1 的维度有助于对网格创建、矩阵和向量运算等进行向量化处理。
假设矩阵 A
代表考试分数,各种行表示不同的班级。您需要计算每个班级的平均分数与各分数的差。使用循环,该运算如下所示:
A = [97 89 84; 95 82 92; 64 80 99;76 77 67;... 88 59 74; 78 66 87; 55 93 85]; mA = mean(A); B = zeros(size(A)); for n = 1:size(A,2) B(:,n) = A(:,n) - mA(n); end
要执行该运算,一种更直接的方式是使用 A - mean(A)
,后者不需要使用循环并且速度明显更快。
devA = A - mean(A)
devA = 18 11 0 16 4 8 -15 2 15 -3 -1 -17 9 -19 -10 -1 -12 3 -24 15 1
即使 A
是一个 7×3 矩阵,mean(A)
是一个 1×3 向量,MATLAB 也会隐式扩展该向量,就好像其大小与矩阵相同一样,并且该运算将作为正常的按元素减法运算来执行。
根据操作数的大小要求,对于每个维度,各数组必须大小相同,或者其中一个为 1。如果满足该要求,则含有大小为 1 的数组的维度将扩展为其他数组中相应维度的大小。有关详细信息,请参阅基本运算的兼容数组大小。
并且,在您处理多维数据时,隐式扩展也有助于向量化操作。假设您要计算包含两个变量 x
和 y
的函数 F
。
F(x,y) = x*exp(-x2 - y2)
要在 x
和 y
向量中的每个点组合处计算该函数,您需要定义网格值。对于此任务,您应避免使用循环来循环访问点组合。在这种情况下,如果一个向量为列,而另一个向量为行,则当这些向量与数组运算符配合使用(例如 x+y
或 x-y
)时,MATLAB 将会自动构造网格。在此示例中,x
是一个 21×1 向量,y
是一个 1×16 向量,因此该运算会通过扩展 x
的第二个维度和 y
的第一个维度来生成一个 21×16 矩阵。
x = (-2:0.2:2)'; % 21-by-1 y = -1.5:0.2:1.5; % 1-by-16 F = x.*exp(-x.^2-y.^2); % 21-by-16
逻辑数组运算
批量处理数组的逻辑扩展是实现比较的向量化并制定决策。MATLAB 比较器接受向量输入并返回向量输出。
例如,假设在从 10,000 个圆锥体中收集数据时,您记录多个负值作为直径的值。您可以使用 >=
运算符确定向量中的哪些值是有效的:
D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14]; D >= 0
ans = 0 1 1 1 0 1 1
Vgood
,其对应的 D
元素是非负的:Vgood = V(D >= 0);
MATLAB 允许您分别使用 all
和 any
函数对整个向量的元素执行逻辑 AND 或 OR 运算。如果 D
的所有值都小于零,则会引发以下警告:
if all(D < 0) warning('All values of diameter are negative.') return end
MATLAB 也可以比较两个大小兼容的向量,并允许您施加更多限制。此代码查找 V 为非负且 D
大于 H
的所有值:
V((V >= 0) & (D > H))
为便于比较,MATLAB 包含特殊值来表示溢出、下溢和未定义的运算符,例如 Inf
和 NaN
。逻辑运算符 isinf
和 isnan
有助于针对这些特殊值执行逻辑测试。例如,往往有助于从计算中排除 NaN
值:
x = [2 -1 0 3 NaN 2 NaN 11 4 Inf]; xvalid = x(~isnan(x))
xvalid = 2 -1 0 3 2 11 4 Inf
注意
Inf == Inf
返回 true;但 NaN == NaN
始终返回 false。
矩阵运算
在使代码向量化时,您通常需要构造一个具有特定大小或结构的矩阵。有多种方式可用来创建均匀矩阵。例如,您可能需要一个包含相等元素的 5×5 矩阵:
A = ones(5,5)*10;
v = 1:5; A = repmat(v,3,1)
A = 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
函数 repmat
可以灵活地根据较小的矩阵或向量来构建矩阵。repmat
通过重复输入矩阵来创建矩阵:
A = repmat(1:3,5,2) B = repmat([1 2; 3 4],2,2)
A = 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 B = 1 2 1 2 3 4 3 4 1 2 1 2 3 4 3 4
排序、设置和计数运算
在许多应用程序中,对向量的某个元素执行的计算取决于同一向量中的其他元素。例如,向量 x 可能表示一个集。如何循环访问一个集而不使用 for
或 while
循环并不明显。如果采用向量化代码,则此过程变得更加清晰而且语法也不再过于冗长。
消除冗余元素
有许多不同的方法可以用来查找向量的冗余元素。一种方法涉及到函数 diff
。在对向量元素进行排序后,当您对该向量使用 diff
函数时,相等的邻接元素会产生零项。因为 diff(x)
生成的向量的元素数比 x
少一个,所以您必须添加一个不等于该集中的任何其他元素的元素。NaN
始终满足此条件。最后,您可以使用逻辑索引来选择该集中的唯一元素:
x = [2 1 2 2 3 1 3 2 1 3]; x = sort(x); difference = diff([x,NaN]); y = x(difference~=0)
y = 1 2 3
unique
函数完成同样的操作:y=unique(x);
unique
函数可能提供超出需要的功能并且降低代码的执行速度。如果您要衡量每个代码节的性能,请使用 tic
和 toc
函数。计算向量中的元素数
您可以计算某元素在向量中的出现次数,而不只是返回 x
的集或子集。该向量排序后,您可以使用 find
函数来确定 diff(x)
中零值的索引并显示元素的值在何处发生了改变。find
函数中后续索引之间的差可指明特定元素的出现次数:
x = [2 1 2 2 3 1 3 2 1 3]; x = sort(x); difference = diff([x,max(x)+1]); count = diff(find([1,difference])) y = x(find(difference))
count = 3 4 3 y = 1 2 3
find
函数不返回 NaN
元素的索引。您可以使用 isnan
和 isinf
函数计算 NaN
和 Inf
值的数目。count_nans = sum(isnan(x(:))); count_infs = sum(isinf(x(:)));