matlab提速matlab提速的问题

下面说说如何对matlab提速的问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ，我会使用两个例子来说明。 【例1】任务描述：根据A的取值使用imshow函数显示矩阵B A = randn(100, 100); B = zeros(size(A)); 【 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 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的多维数组问题”一文中已经提到可以把多维数组转化为一维数组来处理，在这里，也要利用这个思想。需要说明一下，在matlab中，只要循环总次数固定，无论是采用多少个for语句，其运行时间都是不变的，当然了，前提是循环体不变。也就是说，在时间上使用多维数组和一维数组是一样的（不过在空间上我想一维数组要占优），因此，在这里使用一维数组貌似没有什么用处，然而，从下文可以看到，正是因为 “多维变一维”的出现，才令循环的提速得以实现。 首先利用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运行速度慢，不如先改进一下你的算法吧。 之二：（写于2006.09.23） 上一篇综合介绍了如何利用混合编程技术和优化代码技术来提高程序在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(Dsqrt(3*(i-1)^2)) 语句就足够了），但是该函数要消耗巨多的CPU时间，因此只有放弃。 其他具体的细节请参阅我写的m文件（程序有注释）吧，当然了，前提是你有兴趣、耐性和一定的基本功，呵呵。话说回来，用matlab写的程序虽然高度集成，速度快，但是却不方便阅读。同时，也欢迎高手指点。 后记：从我对论坛这一个月的观察来看，基本上大家都是把matlab当作C语言的辅助工具来使用的，出现的帖子更多是询问如何调用matlab程序（例如VC调用com组件）等类似的问题。不过这也无可厚非，谁叫matlab在处理数据方面做得如此出色呢？不过，我还是建议大家认真学习matlab，就算你的初衷是在VC调用matlab的程序，也总不能把代码写得“循环满天飞”，到头来只会抱怨matlab不好，自己的水平也无法提高。 之三：（写于2006.11.19） [原创]提高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那样，有一个团队可以互相讨论互相提高，因此，我的心得难免有错误或理解不当的地方，请大家见谅