日志
提高matlab运行速度和节省空间的心得合集[转]
||
提高matlab运行速度和节省空间的心得合集
2009-01-05 10:51
|
首先推荐使用matlab 2006a版本，该版本优点很多（不过有一个小bug，就是通过GUI自动生成的m文件居然一大堆warning，希望在已经发布了的2006b版本中有改善），其中对于编程人员来说比较突出的一个就是编辑窗口的自动语法检查功能。这可以在一定程度上避免使用没有被定义或赋值的变量，另外，也可以帮助你优化代码，【例1】的【方案3】就是因为我看到matlab编辑窗口的warning而得到的启发。顺便提一下，虽然matlab不像其他语言那样，对变量采用“先定义，后使用”的规则，但是，我的经验是，在使用一个变量之前，最好先对它进行“定义”，这里的“定义”是指为它分配空间，这样不但可以提高运行的速度（这在matlab的帮助中也提到，详见Preallocating Arrays一节），而且还可以减少出错的几率，特别是在循环赋值、且变量大小不固定的时候（对此可参阅这个帖子：[url]http://vib.hit.edu.cn/forum/thread-23732-1-9.html[/url]）。 下面说说如何对matlab提速的问题，我会使用两个例子来说明。 【例1】任务描述：根据A的取值使用imshow函数显示矩阵B A = randn(100, 100); B = zeros(size(A)); 【方案1】 [X,Y] = find(A > 0.6); For i = 1:length(X) B(X(i),Y(i)) = 1; End 【方案2】 B = zeros(size(A)); X = find(A > 0.6); B(X) = 1; 【方案3】： B = zeros(size(A)); X = logical(A > 0.6); B(X) = 1; 事实上，【方案1】到【方案2】的改进在“再谈Matlab的多维数组问题”一文中已经提及过，但是没想到自己在矩阵输入的时候注意到了，但是在矩阵输出的时候就忘记了，看来程序是需要不断修改、优化的，技巧也是需要不断巩固的。至于【方案3】，是得益于matlab的warning提示。 然而，这并不是表示所有类似的地方都可以用logical代替find，当遇到循环次数与X有关时，用find会更有效，对此可以参考【例2】的【方案2】，这是用logical实现不到的。 【例2】任务描述：有一个四维矩阵A（大小为61*73*61*210），其中前三维表示一个包含大脑结构的立方体，最后一维表示大脑中每个点对应的一个长度为210的时间序列。另外有一个三维矩阵Mask（大小是61*73*61），B是二值的，其中1表示该点是前景点（大脑），0表示该点是背景点。任务是对A中属于前景点的时间序列进行EMD处理，从而判断该前景点是否属于激活区。 显然，这个问题要利用循环来完成。在本人写的“再谈Matlab的多维数组问题”一文中已经提到可以把多维数组转化为一维数组来处理，在这里，也要利用这个思想。而且，从下文可以看到，正是因为 “多维变一维”的出现，才令程序得到更进一步的提速。 首先利用reshape函数把四维矩阵A变成二维矩阵B，把三维矩阵Mask变成一维矩阵C： B = reshape(A, 61*73*61, 210); C = Mask(:); t = 1:210; 【方案1】 iTotalVoxel = 61*73*61; for k = 1:iTotalVoxel if C(k) == 1 temp = B(k,: ); imf = emd(t, temp); … end end 【方案2】 D = find(C); iTotalBrainVoxel = length(D); for k = 1: iTotalBrainVoxel temp = B(D(k),: ); imf = emd(t, temp); … end 【方案1】明显是基于C语言的套路，而【方案2】则充分避免了matlab的弱点――循环，经过改进以后（“多维转一维”为此提供了保证，多维的话，恐怕要使用形如A(X(k),Y(k),Z(k),:)的形式了），由于循环次数的降低（大约降低为原来的1/3），故运行时间大致上减少了一半。 由此可见，在matlab中，想加快运行速度，不但要减少循环的层数，而且，还要减少循环的次数。 以下是对【例2】（实现大脑激活区检测）的不同实现的结果比较：（核心工作是对73000个前景点相应的时间序列进行EMD处理，生成10个左右的imf，其中抽取每个imf时平均迭代大概10次左右） 版本1：完全用matlab编写——运行时间大约是200分钟，也就是说，要在matlab做73000*100次循环，最头痛的是迭代时的需要产生样条包络，默认的spline函数耗时相当多； 版本2：用matcom转换成cpp文件后在Borland C++ Builder中运行，完全脱离matlab——由于spline函数耗时太多，因此转换前改用了折线包络，而非样条包络，运行时间为15分钟左右，不过这个结果是对仿真数据，非实际数据而言的，因为我仍未解决matrix.h和matlib.h不能共存的问题，故无法对实际数据进行测试，不过一般来说实际数据比仿真数据运行速度更慢； 版本3：把核心代码（基于样条包络的EMD算法）做成dll文件后在matlab中调用——整个程序，从数据输入到数据输出，只在一个地方使用了一层for循环（就是【例2】中【方案2】的循环），且结合上述优化方案，对于实际数据大概5分钟就能得出结果。 小结：要学好matlab，有效地使用matlab，一定要摆脱C＋＋的思想。能不用循环的地方，尽量不要使用（例如求极值点等这些算法基本上可以不使用循环便可实现），逐渐抛弃C＋＋的逐点运算的思想，多从矩阵的整体（或分块）上考虑。虽然matlab 2006a版对循环已经改进不少，但是，循环的确是造成程序运行速度降低的主要原因。matlab提供的远远不止在调用函数上的方便（例如在C++中编写fft、dwt等函数可能需要几十甚至几百行，而在matlab中只需要一两个语句），运行速度慢或许是没有使用好它，让它发挥出所长所致的。想matlab更高效地为你服务，那就需要不断修改、优化你的代码吧（我的程序编写大概用了一个星期，而修改、优化的时间就用了两个多星期，呵呵）。最后，套用某人的一句话来作结：与其抱怨matlab运行速度慢，不如先改进一下你的算法吧。 [color=Purple][size=2]之二：（写于2006.09.23）[/size][/color] 上一篇综合介绍了如何利用混合编程技术和优化代码技术来提高程序在matlab中的运行速度，这次主要介绍在matlab中如何通过改进算法思想来提高运行速度。 首先介绍一下任务：在图像分割领域，有一种区域增长算法。它的基本思想是将具有相似灰度性质的象素合并成一个较大的区域。已有文章（这是我师兄完成的）介绍可以利用区域增长的思想将fMRI数据中某象素的时间序列与其邻域的具有相似性质的时间序列进行合并，但由于程序丢失了，所以我要重新自己编写。其中fMRI数据是一个时空的四维矩阵，大小为M×N×O×T，前三维代表三维的空间点，最后一维代表每个三维点（voxel，跟二维图像的pixel的意义一样）对应的时间序列。用F = {f_{xyzt}}_{MNOT} 来表示fMRI 数据，其中“_”表示下标，下同。 基于区域增长的激活区检测算法分为区域增长和区域选择两步。在区域增长阶段，图像中的每个象素分别被作为种子象素并依据相似准则进行增长。在区域选择阶段，根据预先定义的筛选规则，某些感兴趣的区域或者任务相关的区域可以被提取出来。 整个算法以伪代码的形式表示，参见附件。 对于附件中的伪代码，有几个地方需要说明一下，如果读者只想理解基本思想可自行跳过这些细节： 1. 集合E 为：所有与区域A_i , i ∈{1,2,...,M × N ×O} 相邻且未被分配的象素，用公式表示为： E = { q ∉ A_i | N(q) ∩ A_i ≠ 0} 其中N(q) 是象素q 的邻居象素。 对于邻居象素，在三维的体图像中，有6－邻域，18－邻域，和26－邻域之分。我用的是26－邻域。 2. 相似性准则的定义： Corr( f_{xyz} ,M_i) > T_1 其中，f_{xyz} 是象素(x, y, z) 对应的时间序列， M_i 是已合并区域 A_i 的平均时间序列。Corr(*,*) 是 Pearson 相关系数，T_1 为预先定义的阈值。该准则matlab提供相应的corrcoef函数实现。 3. “在集合A中找出最大的集”中最大的含义为含有最多的象素个数。 4. {A_j}_J是所有以 A_i 中的元素为种子所生成的集合，即 J ⊂ {1,2,...,M × N ×O}, 其中J定义为 J = { j| v_j ∈ A_i }。 显然，该伪代码可以直接在matlab中实现，其中瓶颈部分是“区域增长”所用到的三层循环，并且最重要的是，在第三层循环中，随着半径的增大（以当前点i为圆心）需要检查的象素会剧增，其增长速度是非线性的，在最坏情况下（每个邻居点都被选中）达到i^2数量级（精确数字是24*i^2+2）！要在matlab中进行这么多次循环其运行速度是很慢的。据师兄的论文记载，运行时间大概是：对于64*64*20*160的仿真数据，需要8分钟；对于61*73*61*210的实际数据，需要260分钟。 下面是我根据算法重新整理后的matlab实现，思想是把第三层循环去掉，从整体去考虑，即只使用两重循环，这样无论三维矩阵的大小如何增大，我的算法都不含有非线性增长速度的因子，顶多是随着i的最大值（即最大半径）的增加而增大循环次数，我想这个循环是无法避免的。改进后的速度还没有具体测试，估计在10分钟内可以做完。改进的基本做法是： 1. 利用matlab的“:”运算符选择邻居象素； 2. 利用matlab的逻辑变量(logical)来进行数组或矩阵的下标引用； 3. 利用matlab的find函数求出数组或矩阵的索引； 4. 利用bwlabeln函数来求得连通区域，对于前后两次半径，如果通过连通区域检测到属于该voxel的区域没有发生变化，则停止增长。 5. 本来想利用matlab的bwdist函数求出其他点跟当前点的距离，以代替求邻居象素的索引（使用D = bwdist(A)后，一个 ind = find(D<sqrt(3*i*2) & D>sqrt(3*(i-1)^2)) 语句就足够了），但是该函数要消耗巨多的CPU时间，因此只有放弃。 其他具体的细节请参阅我写的m文件（程序有注释）吧，当然了，前提是你有兴趣、耐性和一定的基本功，呵呵。话说回来，用matlab写的程序虽然高度集成，速度快，但是却不方便阅读。同时，也欢迎高手指点。 后记：从我对论坛这一个月的观察来看，基本上大家都是把matlab当作C语言的辅助工具来使用的，出现的帖子更多是询问如何调用matlab程序（例如VC调用com组件）等类似的问题。不过这也无可厚非，谁叫matlab在处理数据方面做得如此出色呢？不过，我还是建议大家认真学习matlab，就算你的初衷是在VC调用matlab的程序，也总不能把代码写得“循环满天飞”，到头来只会抱怨matlab不好，自己的水平也无法提高。 [color=Purple][size=2]之三：（写于2006.11.19）[/size][/color] [原创]提高matlab运行速度和节省空间的一点心得（之三） 首先说说Matlab与其他语言的差异：例如对于C或者C++来说，只要算法的思想不变、采用的数据结构相同，不同人写出来的语句在效率上一般不会产生太大的差别。所以，对于C来说，程序的好坏一般由算法来决定。但是，在matlab中，同样的算法、同样的结构、同样的流程，如果采用的语句不一样，在效率上就会大大不同。所以，我认为，使用matlab比使用其他语言更加困难，也显得matlab更难以掌握。另外，由于matlab在存储管理上的不便，使得在同时提高时空两域的效率变得更加困难，特别是在空间上（因为在时间上matlab提供了profiler这个非常有用的工具，但是在空间上就没有）。当需要处理大量的数据时，精简时空两域的程序语句就尤为重要了。 空间上： 1. 建议使用A = logical(sparse(m,n))，不建议使用 A = sparse(false(m,n))，两者结果一样，但是后者生成m×n的临时矩阵，浪费空间，且当m、n很大时，后者不一定能申请成功； 2. 使用sparse几点注意： a) 只能用在二维以下的矩阵上； b) 由于matlab按照“先行后列”的方式读取数据（即先把第一列所有行读取完以后再读取第二列的各行），因此定义稀疏矩阵时，最好“行数>列数”，这样有利于寻址和空间的节省（自己试试a=sparse(10,5); whos a和b= sparse(5,10);whos b就知道了）； c) 对大型矩阵用sparse非常有效（不但节省空间，而且加快速度，强烈推荐！这在动态申请数组空间的时候尤其方便，当然了，数组不是太大的时候也可以使用eval即字符串的方法），但对小型矩阵使用反而增加存储量（自己试试a=false(5,1); whos a和b=logical(sparse(5,1));whos b就知道了），相信这是由稀疏矩阵需要存储额外的信息引起的。 3. 尽量按照精度来选择数据的类型，例如，如果只需用到0－255之间的整数，则定义矩阵为uint8型就ok了，定义方式：A = zeros(10,10,’uint8’)；可以用intmin(‘uint8’)和intmax(‘uint8’)返回该种类型的最值。 时间上： 1. 在for循环中，清零操作用赋值语句 A = B，其中B是在for循环外的一个同A大小一样的全0阵，不要使用A(:) = 0；但这样会大大影响后面的逐点运算速度。这个问题要请教高手，那就是“个别语句的改动会引起其他语句的执行速度”。例如分别运行3万多次下面代码，但执行时间有较大差别： iLen = length(find(alRegion)); % 0.58 s if iLen >= iThreshold_2 % 0.05 s … end 和 iLen = sum(alRegion); % 0.37 s if iLen >= iThreshold_2 % 0.40 s … end 2. Find函数较慢，可用logical函数代替，但是，当需要取得满足条件的下标时，就无法使用logical函数，这已经在我之前的帖子中提及过。不过，大家有没有想过，连find函数都可以进行优化的，方法是用“基本矩阵进行显式逻辑引用”来代替find。例如，假设矩阵A是一个61*73*20的三维逻辑矩阵，如果用下面的语句循环3万多次，需要的时间是13 s ： B = find(A == true); 如果采用下面的方法，则只需要不到0.7 s： 首先定义一个索引矩阵： a3iCubicIdx = uint32(1:iTotalVoxel); % uint32可以根据需要调整，这里省略了条件判断 a3iCubicIdx = reshape(a3iCubicIdx, [iVolLen, iVolWdh, iVolHgh]); 然后在循环中写以下代码： a3iTemp = a3iCubicIdx(iXmin:iXMax,iYMin:iYMax,iZMin:iZMax); B = aiTemp(A(iXmin:iXMax,iYMin:iYMax,iZMin:iZMax)); 当然了，改进的前提是知道矩阵A的非零元（即值为true的元素）大致的分布，也就是能够求出iXmin:iXMax,iYMin:iYMax,iZMin:iZMax这个范围。现在终于明白并体会到cwit所说的“连num2str都优化过”的含义了。 3. 不断优化代码，例如corrcoef函数，matlab自带的corrcoef函数求整个矩阵所有列的相关系数，因为我只需要求出某一列跟其他各列的相关系数，所以参照corrcoef函数自己写了一个，不但把速度提了上去，而且还发现了：repmat(5，100，1)的速度并不比ones(100,1)*5 快，另外，别小看一个小矩阵的转置操作，当循环次数很大的时候，有没有转置就相差很远了。 4. 使用逻辑运算符&、| 时，两个操作对象最好是logical类型，否则速度会减慢。 5. 二维矩阵转置操作可以用以下三种方法进行，三者的效率基本一样（时空），如果遇到三维以上的矩阵要转置，用permute命令较为方便： a) A = A’; b) A = permute(A,[2,1]); c) A = shiftdim(A,1); 6. “使用eval方式动态存储多个一维数组”比“使用二维数组动态存储多个一维数组”要快，即：eval(['A_', num2str(k),' = B;']);比 A(k,:) = B; 快，其中B是一个一维数组，k表示循环次数。注：并非所有B都进行存储，只存储满足某个条件的B，另外，对预申请空间A不成功，这是对结论的补充说明。值得注意的是，如果对B是一个稀疏的一维数组，则eval方式的优势荡然无存，当k增大时反而增加系统开销。 7. 当矩阵很大时，利用A(:,k+1:end)=[];去掉多余元素操作时会减慢程序的运行，因此，如果后续处理中没有用到这些多余元素，则没有必要使用这个语句，即不管就是了。 8. 当需要对很大的一个矩阵进行操作时，可以考虑使用循环来完成。例如corrcoef函数，如果处理的对象矩阵A是100*180时（即对100个列向量求它们两两之间的相关程度，假设需要的只是前面99个与第100个向量的相关系数，其他不需要用到），直接用corrcoef(A)会比较慢，这时候可以考虑把矩阵A分为5个部分，每个子块与第100个向量进行相关，这样速度会更快。 9. 局部比较、赋值比全局比较、赋值要快（呵呵，这是废话），假设A、B都是三维逻辑矩阵，如果只想对某个局部（例如X_1:X_2,Y_1:Y_2,Z_1:Z_2这个立方体）进行比较和赋值，则推荐使用B(X_1:X_2,Y_1:Y_2,Z_1:Z_2) = B(X_1:X_2,Y_1:Y_2,Z_1:Z_2) & ~A(X_1:X_2,Y_1:Y_2,Z_1:Z_2)，这比B(A) = false或B = B&~A速度上都要快不少。 后记： 1. 此心得原来打算在cwit给了我回复，指正文中错误之后再整理并发布的，不过cwit很忙，所以不知道要等到什么时候，因此先献丑了，有错误大家也帮忙指正一下。（注：这个版本删除了“循环总次数固定无论使用多少个for循环速度不变”这句，因为考虑到矩阵维数会影响速度） 2. 上述技巧都是出自我在处理大型矩阵时候自己总结出来的个人心得，转载时请注明。 3. 实验室中由于只有我一个人比较关心matlab的运行速度和存储空间问题（尤其是速度问题），所以不像cwit那样，有一个团队可以互相讨论互相提高，因此，我的心得难免有错误或理解不当的地方，请大家见谅。 [[i] 本帖最后由 eight 于 2007-1-16 13:45 编辑 [/i]] eight发在振动上的原创贴，一共三贴，前后时间贯穿三个月，有值得借鉴的东西。应我请求，合在一起贴在simwe，是总结了很长时间的心得，欢迎各位接续讨论。我把在振动上我的帖子也转过来，其中观点未必正确，同样仅供参考： [color=brown][b]1.你计算时间是用profiler还是tic/toc？两个你想比较速度的程序中，第一个运行完毕后你是否关闭MATLAB还是直接运行第二个(如果单单用clear命令我认为上述结论就未必合理，clear的方式对内存的清理不是很彻底，和MATLAB挂钩的比如计算engine是两个程序相同的，clear并不能清除这些东西)？ 2.eval命令我个人不赞成大量使用，string翻译成命令再执行，速度下降相当明显，文章似乎未曾提到这一点。 3.repmat命令感觉如果单纯建立全1阵时肯定不会比ones快，但是功能高于ones，毕竟是复制matrix的命令。 4.部分内容已经有些讨论，比如find命令肯定没有矢量化的索引方式更快等。 5.补充：三维矩阵没事儿尽量不要用，改cell或者struct代替，可以去测试，系统资源的占用远远大于后两者。[/b][/color] [[i] 本帖最后由 bainhome 于 2006-11-25 22:56 编辑 [/i]] [quote]原帖由 [i]bainhome[/i] 于 2006-11-25 22:54 发表 eight发在振动上的原创贴，一共三贴，前后时间贯穿三个月，应有很多值得借鉴的东西。应我请求，又合在一起贴在simwe，应该是总结了很长时间的，欢迎各位接续讨论。我把在振动上我的帖子同样也转过来，其中观点未必 ... [/quote] 1. 计算时间用的是profiler，时间统计是把gui（我所指的gui是指file--new--gui那个方式）中的核心程序单独放在一个m文件中，然后运行一次完毕后，只修改那个地方后再运行统计的，没有关闭matlab 2. eval命令在矩阵size小时比较快捷的，按照我的运行结果，100次左右时比二维矩阵快，当然了，你说得无错，size大时系统开销很大，我在文中也有提及。事实上，eval命令是我在使用稀疏矩阵前采用的方法，后来发现稀疏矩阵后，就没有用eval了 3. 三维矩阵有些时候无法避免，例如我要处理的是人脑坐标，大小是61*73*61或者是79*95*68，而且，不能把它变成一维，因为要用到三维空间信息 先支持一下，稍后再参与讨论。 你说得有些情况，在开发工程软件时，条件会变，结论也会变。 赞一个！俺也贴两条，可能众所周知:D 1. 多重循环中把循环次数多的放内层有利于提高计算速度。 交换内外循环次数试试下述代码： tic; for n=1:10 for m=1:1000000 sum = m+n; end end toc 2. 用结构数组（Struct Array）比用数组结构（Array Struct）节省内存。 S1.R(1:100,1:50) = 1; S1.G(1:100,1:50) = 2; S1.G(1:100,1:50) = 3; for m=1:50 for n = 1:100 S2(n,m).R = 1; S2(n,m).G = 2; S3(n,m).B = 3; end end whos S1 S2 S1和S2存储的数据一样，但S1用得内存要比S2少得多。 *************************************************************** ** 建议大家继续跟贴，把自己在优化速度和内存方面的心得拿出来分享！** *************************************************************** [[i] 本帖最后由 junziyang 于 2006-11-29 09:52 编辑 [/i]] 我也帖一个 matlab列优先，所以外循环列，内循环行，要快，提高了在cache的命中率 good roger! that's cute 在MATLAB运用方面我还是个新手,暂没有相关方面的心得拿出来分享 以后会常来这里学习 谢谢分享原创经验! 要是加上运行时各方案的具体耗时就更加生动了。 7. 当矩阵很大时，利用A(:,k+1:end)=[];去掉多余元素操作时会减慢程序的运行，因此，如果后续处理中没有用到这些多余元素，则没有必要使用这个语句，即不管就是了。 关于这句 不是很明白? 去掉反而减慢速度.而且很后面用倒有什么关系 ? :victory: 好方法！ |
