MATLAB如何尽量避免for循环？

我们知道MATLAB擅长矩阵计算，但对于跑for循环非常低效，因此在内存足够的情况下应尽量写成矩阵或者向量化操作的形式，善用更好的数据结构、算法，以及matlab自带的函数特性，以尽可能避免for循环降低运行速度。下面是我学到的一些小tips，并结合运算示例进行讲解。

方法1：优化循环嵌套及内部运算

例如：

1. 将计算量小、循环次数多的放在里面，计算量大循环次数少的放在外面做大循环；
2. 将表达式尽可能向量化，计算好后再放入循环内，减少循环內部运算，避免循环中调用计算量很大的函数；

方法2：向量化处理

1. 初始化变量，即预分配内存以获得更高运算速度，避免变量在循环迭代中改变。这一类的函数有zeros, ones, cell, struct, repemat等。注意选择合适的数据类型。如果能满足需求，分配内存时尽量用更低精度数据类型例如logical, uint8, int8等以节省内存和加快运算（注意若涉及浮点数运算还是得用float或者double）。例如预定义加载到DMD上的binary pattern，可以直接设置成logical类型：

DMDpattern = zeros(1920, 1080, &＃39;logical&＃39;);


2. 将循环语句改为向量化形式，关键是构造索引数组，（在不同维度上）同时索引矩阵多个元素。

• 例如要画出 $y=sin(x)$图像，MATLAB许多数学函数都支持向量化输入，直接写成如下形式，而不要像C语言那样依次循环赋值给因变量。

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


• 例如要对矩阵的偶数列进行上下翻转，奇数列进行上下移位，可以不用循环，直接构造索引数组，example code如下：

for i=2:2:P_x %矩阵翻转，偶数列x轴扫描时回来的数据reverseOR_MAP(:,i)=flipud(OR_MAP(:,i));
end
for i=1:2:P_x %由于电机奇数列扫瞄过去时可能错位，需往上平移，否则图像抖动shiftY = 4;for j = max(shiftY+1,1):1:min(P_y, P_y+shiftY)OR_MAP(j-shiftY,i) = OR_MAP(j,i);end
end


改为：

%矩阵翻转，偶数列x轴扫描时回来的数据reverse
evenIdx = 2:2:P_x;
OR_MAP2(:,evenIdx) = flipud(OR_MAP(:,evenIdx));
%由于电机奇数列扫瞄过去时可能错位，需纵向平移，否则图像抖动
oddIdx = 1:2:P_x;
shiftY = 4;
OR_MAP((max(shiftY+1,1):min(P_y, P_y+shiftY))-shiftY, oddIdx) = ...OR_MAP(max(shiftY+1,1):min(P_y, P_x+shifty), oddIdx);


方法3：善用MATLAB向量化函数特性

向量化(vectorizing) 操作这个过程将向量扩展为更大的索引矩阵，以尽可能替代for循环，本质是用空间换时间。MATLAB中有一些支持vectorizing技术的常用函数如：All, diff, ipermute, permute, prod, reshape, squeeze, any, find, sort, sum, repemat, repelem等， 更高级的还有 kron, shiftdim, cumsum, ndgrid, sub2ind in2sub, bsffun, arrayfun 等函数，充分利用这些函数特性非常重要，下面我们举一些应用实例：

1. 矩阵扩充：使用repmat函数复制矩阵，repelem函数复制矩阵元素，kron函数计算矩阵Kronecker积
2. 用 find 函数寻找满足要求的矩阵元素索引，例如需要找到矩阵ORdata每列的最大值所在的行序号，其中矩阵某列最大值可能有多个而我们只需要找到第一个最大值 ，example code如下：

[ORmaxRow, ORmaxCol] = find(ORdata==max(ORdata));
[~, ia, ~] = unique(ORmaxCol);
ORTmax = ORmaxRow(ia);


3. sub2idx函数用于将多维数组索引转化为特定维度矩阵的线性下标，以便赋值

例如：有一个矩阵，上面有索引和需要赋予的数值，需要把对应的数值赋予到三维数组中，当时由于数据量巨大（数万个点），有没有可以不用for循环的方法进行赋值？
举个例子，比如: A=[1 1 2 42;2 2 3 997;3 4 5 10000;4 4 7 899];
已存在一个大小为mxnxq的三维数组B，以A矩阵的前三列为索引，把第四列的数值赋予到三维矩阵中，在这个例子中，结果就是
B(1,1,2)=42 B(2,2,3)=997 B(3,4,5)=10000 B(4,4,7)=899

如何不用for循环做到这个效果？

clear;
A=[1 1 2 42;2 2 3 997;3 4 5 10000;4 4 7 899];
B=zeros(4,4,7);B(sub2ind(size(B),A(:,1),A(:,2),A(:,3)))=A(:,4);


4. 使用 ind2sub将线性下标转换为特定维度矩阵的坐标数组

例如：在三维坐标系下（X，Y，Z），得到一个三维数组E（20,20,20），也就是X，Y，Z各维度上有20个点，需要把它化成一维数组，并且得到相应的索引（坐标）数组。

直接E(:)或者用reshape，得到的只是形状变成（800,1）的原数组。我们可以用 ind2sub提取X, Y, Z坐标数组。ind2sub用法如下：

IND = [3 4;5 6];
s = [2,2,2];
[I,J,K] = ind2sub(s,IND)


5. 构造匿名函数结合 arrayfun 从大型矩阵中同时提取多个子矩阵

例如：从二维矩阵cam中提取15个不同的100×100矩阵块。这可以使用for循环轻松完成，但如果不用循环的话，有什么优雅的解决方案吗？

错误示例：可以构造包含矩阵块左上角索引向量x和y，如下代码将只产生100×100矩阵而不是15x100x100，因为x(1:end):(x(1:end)+99)仍然还是一维数组。

result = cam(x(1:end):(x(1:end)+99), y(1:end):(y(1:end)+99));


正确做法：我们可以构造匿名函数结合 arrayfun 和 cat，代码如下：

indexFcn = @(r,c) cam1(r:(r+99),c:(c+99));
result = arrayfun(indexFcn,x,y,&＃39;UniformOutput&＃39;,false);
result = cat(3,result{:});


说明：

第一行创建一个匿名函数，这是一个简单的单行函数，可以在运行时创建，而无需将其放在m文件中。该函数定义了两个输入r和c，用于从cam中提取100×100的子矩阵.变量indexFcn存储用于调用该函数的function handle。

第二行调用 arrayfun,它将函数应用于数组的每个元素. arrayfun 循环遍历x和y中的每个条目,将值传递给indexFcn.输出存储在结果中，这是一个15元素的元胞阵列,其中每个元胞包含一个100x100的矩阵.

第三行使用cat函数将包含15个100×100矩阵的元胞阵列连接成100×100×15的矩阵.

6. 利用permute 函数重排矩阵维数，shiftdim 函数偏移矩阵维数，实现巧妙的矩阵运算

例如：A=[a b c d e],
A中都是列向量，想依次计算aa’,bb’,cc’,dd’,ee’,储存到一个三维矩阵B(:,:,5)，不用for循环一个一个计算出来，该怎么做？实际计算的矩阵有13列

permute(A, [1 3 2]) .* shiftdim(A, -1)


说明：

• permute(A, [1 3 2])将2D矩阵A(4x4)转换为每个切片是原矩阵每列的3D数组(4x1x4)
• shiftdim(A, -1)将2D矩阵A(4x4)转换为每个切片是原矩阵每列转置的3D数组(1x4x4)，
• 事实上上述两步也可以通过reshae实现。
• 这样两者做对应元素的数乘得到(4x4x4)的3D数组，其中每个切片就和(4x1)矩阵乘以(1x4)矩阵结果相同，这里.*利用到了MATLAB中的broadcast功能自动拓展。
• .*的自动拓展功能也可以用bsxfun函数替代，C = bsxfun(fun,A,B)通过函数句柄fun能对数组A和B执行pixel-wise的二元操作，如下：

bsxfun(@times, permute(A, [1 3 2]), shiftdim(A, -1))


7. 原本循环实现序列矩阵与向量对应元素点乘，可将序列矩阵reshape为二维，向量变成对角矩阵后进行矩阵相乘

例如原本用for循环，实现3维矩阵T1的每个slice都与1维向量SLM_mask对应元素点乘后求和

E_exc = zeros(size, size);for j = 1:numE_exc = E_exc + T1(:,:,j)*SLM_mask(j);end


• 可以利用reshape函数将序列矩阵压缩为2维，1维向量SLM_mask通过diag转为对角矩阵，矩阵相乘使得二维化的T1每个列向量都乘以SLM_mask对应元素，矩阵列向求和再reshape回来

E_exc = sum(reshape(T1, [], num) * diag(SLM_mask(:)), 2);E_exc = reshape(E_exc, size, size);


参考

[1] 解决Matlab当中for循环运行慢的问题
[2] matlab如何根据索引不用for循环给三维数组赋值？
[3] MATLAB: Extract multiple parts of a matrix without using loops
[4] Matlab中的for循环优化
[5] matlab 如何不用循环完成如下计算?