MATLAB 代码为什么运行得慢从向量化计算到程序性能优化目录MATLAB 代码为什么运行得慢从向量化计算到程序性能优化一、从一个简单的计算问题说起二、MATLAB 为什么适合向量化计算数组运算与矩阵运算向量化不只是使用点运算符三、循环、预分配与向量化三个需要区分的问题数组预分配解决了什么问题向量化改善了什么现代 MATLAB 的循环还慢吗四、隐式扩展MATLAB 数组计算的一个重要特性五、为什么向量化不一定更快不必为了消除循环而滥用 arrayfun六、如何判断一段 MATLAB 代码是否真的需要优化使用 Profiler 查找性能瓶颈七、从向量化回到 MATLAB 编程本身参考资料MATLAB 是一门以矩阵和数组计算为特色的编程语言。对于数值计算、数据处理和科学研究中的许多问题我们往往只需要几行代码就能完成原本较为复杂的运算。然而代码能够正确运行并不意味着它具有良好的执行效率。同一个数学问题采用不同的 MATLAB 实现方式运行时间和内存消耗可能存在明显差异。在讨论 MATLAB 性能优化时经常能看到一条建议尽量使用向量化计算减少for循环。这一建议有其历史背景也确实适用于许多常见的计算任务。不过随着 MATLAB 执行引擎不断改进循环语句的性能已经发生了很大变化。一些原本被认为效率低下的写法在现代 MATLAB 中未必仍然存在明显的性能问题。与此同时过度追求向量化还可能引入额外的内存开销。因此理解 MATLAB 的性能特点需要同时考虑数组运算机制、程序执行方式和数据的存储成本。我们可以从一个简单的数值计算问题开始。一、从一个简单的计算问题说起假设需要在区间 ([0,10]) 上选取大量采样点计算函数的数值。对于熟悉一般编程语言的读者最直接的思路是逐个处理采样点N 1000000; x linspace(0, 10, N); for k 1:N y(k) x(k)^2 * exp(-x(k)); end这段程序在计算逻辑上没有问题。循环依次取出数组x中的元素代入函数表达式并将结果存储在y中。但如果仔细观察会发现这里有两个值得讨论的地方。其一程序没有预先为y分配存储空间。随着循环推进结果数组会逐渐增长。虽然 MATLAB 能够自动管理数组内存但这种写法可能带来额外的内存管理开销。其二所有采样点之间相互独立。计算某个点的函数值并不依赖其他采样点的结果。对于这样的任务没有必要在程序层面逐个组织计算。MATLAB 允许直接将整个数组代入数学表达式N 1000000; x linspace(0, 10, N); y x.^2 .* exp(-x);两段代码完成相同的数学任务第二种写法却更加接近原始函数表达式。这正是 MATLAB 向量化计算的一个典型例子。需要注意代码更短并不自动意味着运行速度一定更快。我们还需要区分数组预分配与向量化的作用并通过实际测量判断不同实现方式的效率。二、MATLAB 为什么适合向量化计算MATLAB 的名称源于Matrix Laboratory矩阵计算从早期就是这门语言的重要特色。与需要显式逐个访问元素的编程方式相比MATLAB 提供了丰富的数组运算符和内置数值函数可以直接对成组的数据进行操作。例如对一个数值向量进行平方、求和或者指数运算x 1:10; a x.^2; b sum(x); c exp(x);这些操作都不需要用户显式编写循环。MathWorks 将向量化Vectorization描述为将基于循环、面向标量的程序改写为使用矩阵和向量运算的程序。其优点主要体现在数学表达方式、代码简洁性以及许多计算场景中的运行效率。数组运算与矩阵运算理解 MATLAB 向量化时必须先区分数组运算和矩阵运算。考虑两个矩阵A [1 2; 3 4]; B [5 6; 7 8];表达式C A .* B;计算对应位置元素的乘积即而C A * B;表示矩阵乘法二者具有不同的数学含义。类似地表示对数组中的每个元素求平方而表示矩阵平方即计算。在处理实验数据、采样序列或离散函数时逐元素运算十分常见。初学者如果忽略这些运算符之间的区别很容易得到维度错误甚至计算出与预期完全不同的结果。向量化不只是使用点运算符除了.*、./和.^MATLAB 还提供了大量能够直接处理数组的函数例如sum、mean、diff、cumsum和sort。假设需要计算一个序列在不同位置的累积和可以使用循环x 1:10000; s zeros(size(x)); for k 1:numel(x) s(k) sum(x(1:k)); end这里每次循环都会重新计算前面若干个元素的和大量工作被重复执行。如果直接使用 MATLAB 的累积求和函数s cumsum(x);代码不仅更加简单计算方法也发生了实质性变化。前一种实现重复计算多个前缀和运算量随序列长度近似呈二次增长后一种方法通过累积计算可以避免这些重复工作。这说明MATLAB 的性能优化有时涉及算法层面的改进。选择合适的内置函数往往比单纯将循环改写成数组表达式更有意义。三、循环、预分配与向量化三个需要区分的问题回到前面的函数计算我们可以构造三种实现方式。第一种实现采用动态扩容function y calcGrow(x) y []; for k 1:numel(x) y(k) x(k)^2 * exp(-x(k)); end end第二种实现保留循环但提前分配结果数组function y calcPrealloc(x) y zeros(size(x)); for k 1:numel(x) y(k) x(k)^2 * exp(-x(k)); end end第三种实现直接使用向量化表达式function y calcVector(x) y x.^2 .* exp(-x); end这三种实现有助于我们区分两个不同的优化来源。数组预分配解决了什么问题MATLAB 的数值数组通常需要连续的内存空间。当数组逐渐增长时执行环境可能需要寻找新的存储区域并移动已有数据。对于一个预先知道结果大小的计算任务可以使用zeros、ones或其他初始化函数提前分配所需空间从而减少数组扩容带来的开销。例如y zeros(1, N);预分配的作用是避免结果数组在循环中不断调整大小。它并没有改变循环的数学计算过程。有必要补充的是现代 MATLAB 对数组增长已经进行了多次优化。早期版本中关于动态扩容的极端性能结论不能直接套用到新版本。MathWorks 的技术文章也指出R2011a 之后自动数组增长的处理方式已经改进。因此预分配仍然是合理的编程习惯但实际收益应当通过测试确定。预分配也不是越多越好。如果后续计算会直接生成一个新的完整数组那么提前创建一个马上被覆盖的零矩阵只会产生不必要的工作。向量化改善了什么对于前面的函数计算各个采样点之间不存在数据依赖因此可以将逐元素操作交给 MATLAB 的数组运算机制处理。这样做减少了用户代码中的显式逐元素循环并能够利用 MATLAB 对数组运算的专门优化。对于许多常见的数学函数和批量数据操作向量化通常具有良好的性能。不过向量化并不意味着计算机在底层不再使用循环。数组中的所有元素仍然需要按照相应的数学规则进行处理只是这些处理过程由 MATLAB 执行引擎及其调用的底层计算例程负责组织。同样向量化也不等于并行计算。某些内置数值运算可能利用多线程计算资源但不能据此推断所有数组表达式都会自动并行执行。现代 MATLAB 的循环还慢吗关于 MATLAB 循环效率的讨论需要考虑软件版本。MATLAB 很早就具有即时编译Just-In-Time CompilationJIT机制。R2015b 发布时MathWorks 对 MATLAB 执行引擎进行了重要重构进一步改进了代码执行和即时编译能力。这一变化意味着现代 MATLAB 已经能够对许多循环、函数调用和逐元素计算进行有效优化。因此“MATLAB 中的循环一定很慢”并不是一个可靠的普遍结论。如果某个算法具有明确的递推关系例如那么后一个状态依赖前一个状态。在没有找到合适的等价数学变换之前使用循环通常是直接且清楚的实现方式。对于这类问题强行消除循环可能导致程序更加复杂却未必带来速度优势。四、隐式扩展MATLAB 数组计算的一个重要特性除了直接对同尺寸数组进行运算MATLAB 还支持维度兼容数组之间的隐式扩展Implicit Expansion。这一机制在 R2016b 中被引入到许多基本数组运算中使得原本需要显式构造重复数据的操作可以直接完成。例如假设矩阵A有 1000 行、5 列每一列代表一个变量需要分别减去该列的均值A rand(1000, 5); mu mean(A, 1); B A - mu;这里A是一个1000×5矩阵而mu是一个1×5行向量。在减法运算中MATLAB 根据数组维度的兼容规则进行隐式扩展使每一行分别减去对应的列均值。对于两个数组只要对应维度的长度相同或者其中一个长度为 1通常就可以进行支持隐式扩展的逐元素运算。在较早的 MATLAB 版本中类似操作可以使用bsxfunB bsxfun(minus, A, mu);相比之下隐式扩展使表达式更加简洁也更接近矩阵数据处理的数学描述。但这一特性同样需要谨慎使用。因为两个数组即使尺寸不同也可能满足隐式扩展条件程序有时会生成比预期大得多的结果。例如x (1:1000); y 1:1000; D x - y;其中x是列向量y是行向量。两者相减会得到一个1000×1000的矩阵而非长度为 1000 的向量。这在构造差值矩阵、二维函数采样和距离计算时非常方便但如果只是因为不小心改变了数组方向就可能产生意外的内存消耗。五、为什么向量化不一定更快向量化的一个重要限制来自内存使用。考虑两个长度均为 (N) 的向量假设需要计算利用隐式扩展可以直接写成s sum(abs(x - y), 2);这段代码很简洁但表达式x-y对应一个的差值矩阵。如果采用double类型存储那么仅一个完整的差值矩阵理论上就需要也就是约 800 MB 的数据空间。计算过程中还可能出现其他临时数组因此实际内存需求可能进一步增加。MATLAB 是否能够省略或复用某些临时存储则与具体表达式及执行优化有关。如果只关心每一行的累加结果而无需保留整个差值矩阵就可以考虑分块计算。下面的实现每次处理部分行N numel(x); s zeros(N, 1); blockSize 200; for first 1:blockSize:N idx first:min(first blockSize - 1, N); D abs(x(idx) - y); s(idx) sum(D, 2); end假设x和y分别为长度为 (N) 的列向量与行向量那么每次循环只需要处理不超过blockSize×N的差值数据。对于blockSize200、N10000的情况单个双精度差值数组的理论存储量约为 16 MB明显小于一次性构造完整矩阵所需的 800 MB。这种写法仍然使用for循环但每个循环内部利用向量化完成一批数据的运算。分块方法的主要优势是可以控制中间数组的规模并非一定能够缩短执行时间。循环次数、内存访问、临时数组分配以及底层运算效率都会影响最终结果。对于较小的数据规模一次性向量化反而可能更简单、更快。由此可以看出循环与向量化并不矛盾。实际程序完全可以同时使用两者通过循环管理计算任务通过向量化处理适当大小的数据块。不必为了消除循环而滥用 arrayfunMATLAB 还提供了arrayfun、cellfun等函数可以将指定函数应用到数组或元胞数组的元素上。这类函数有时能够使程序结构更加紧凑但使用arrayfun并不意味着计算自动获得了更高的执行效率。对于普通 CPU 数组arrayfun仍然需要逐元素应用指定函数具体性能取决于函数调用开销和计算内容。GPU 数组上的arrayfun具有不同的实现机制应当另行讨论。对于简单的数值操作直接采用内置数组运算通常更自然对于复杂逻辑也没有必要仅仅为了避免出现for关键字而改变实现方式。六、如何判断一段 MATLAB 代码是否真的需要优化性能优化不应只依赖经验判断。最直接的方法是在相同测试条件下测量不同实现的实际运行时间。MATLAB 提供了tic、toc和timeit等计时工具。tic与toc使用方便适合测量一段程序运行所需的时间tic y x.^2 .* exp(-x); toc但对于执行时间很短的代码单次计时容易受到系统调度、后台任务及其他因素影响。相比之下timeit会重复调用指定函数并根据多次测量结果给出典型执行时间。对于需要比较不同函数实现的场景通常更适合使用它。前面的三个计算函数可以按照下面的方式测试x linspace(0, 10, 1000000); f1 () calcGrow(x); f2 () calcPrealloc(x); f3 () calcVector(x); t1 timeit(f1); t2 timeit(f2); t3 timeit(f3); fprintf(动态扩容循环: %.6f s\n, t1); fprintf(预分配循环: %.6f s\n, t2); fprintf(向量化计算: %.6f s\n, t3);其中三个函数需要分别定义并保证它们执行相同的计算任务。为了确认优化没有改变计算结果还应当比较不同实现的输出y1 calcGrow(x); y2 calcPrealloc(x); y3 calcVector(x); err1 max(abs(y1 - y3)); err2 max(abs(y2 - y3)); disp([err1, err2]);对于浮点数计算比较结果时应当考虑适当的误差容限而不能不加区分地要求所有算法实现都具有完全相同的二进制结果。在正式性能测试中至少需要记录 MATLAB 版本、计算机处理器、内存配置、输入数据规模及计时方式。如果比较对象涉及不同算法还需要检查其数学结果和数值精度是否满足相同要求。此外只测试一种输入规模通常不足以支持一般性结论。建议分别考察较小、中等和较大的数组观察执行时间随问题规模变化的情况。使用 Profiler 查找性能瓶颈对于包含多个函数和复杂流程的程序逐段使用tic和toc往往不够方便。MATLAB 提供了 Profiler用于分析函数调用和代码执行时间。可以通过以下命令开启和查看性能分析profile on % 在此处运行需要分析的程序 profile off profile viewerProfiler 可以帮助定位耗时较多的函数以及需要进一步检查的代码位置。不过性能分析器本身也需要消耗计算资源。启用 Profiler 后测得的运行时间可能高于程序正常执行时的时间。因此Profiler 更适合用于定位瓶颈关键实现之间的精确计时比较则可以另行使用timeit。这一点非常重要因为实际工程程序中的主要耗时位置未必与我们最初的直觉一致。例如一个程序可能包含大量循环但真正的耗时集中在文件读写、重复的数据类型转换、矩阵分解或某几个高开销函数调用中。此时单纯修改循环写法可能几乎无法改善整体性能。优化工作应当优先针对实际瓶颈而不必让每一行代码都追求理论上的最快执行速度。七、从向量化回到 MATLAB 编程本身向量化是 MATLAB 语言设计中很有代表性的特点。它使数组运算能够采用接近数学表达式的形式编写也为许多数值计算任务提供了简洁而高效的实现方式。但随着 MATLAB 执行引擎的发展将所有循环都视为低效代码已经没有充分依据。对于相互独立的批量计算可以优先考虑数组运算和内置函数对于具有递推关系的算法清晰的循环结构往往更加合适对于数据规模较大的任务还需要关注中间数组和内存分配必要时采用分块计算。实际编写 MATLAB 程序时我更倾向于先选择逻辑清楚、结果可靠的实现方式然后利用计时和性能分析工具识别真正的瓶颈。如果向量化能够让代码更加简洁同时改善运行效率自然值得采用如果某个循环已经足够高效也没有必要为了形式上的向量化而增加程序复杂度。毕竟程序优化的目的始终是以合理的计算资源完成所需任务而不是尽可能减少源代码中的循环语句。参考资料以下资料均来自 MathWorks 官方文档或官方技术博客。1. MathWorks. Vectorization. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html2. MathWorks. Preallocation. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html3. MathWorks. Techniques to Improve Performance. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html4. MathWorks. Compatible Array Sizes for Basic Operations. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/compatible-array-sizes-for-basic-operations.html5. MathWorks. Strategies for Efficient Use of Memory. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/strategies-for-efficient-use-of-memory.html6. MathWorks. Measure the Performance of Your Code. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/measure-performance-of-your-program.html7. MathWorks. Profile Your Code to Improve Performance. MATLAB Documentation.https://www.mathworks.com/help/matlab/matlab_prog/profiling-for-improving-performance.html8. Bergstein, D. Run Code Faster With the New MATLAB Execution Engine. Loren on the Art of MATLAB, February 12, 2016.https://blogs.mathworks.com/loren/2016/02/12/run-code-faster-with-the-new-matlab-execution-engine/9. Eddins, S. MATLAB Arithmetic Expands in R2016b. Loren on the Art of MATLAB, October 24, 2016.https://blogs.mathworks.com/loren/2016/10/24/matlab-arithmetic-expands-in-r2016b/
