记录一些常用的 MATLAB 编程技巧/规范，目标是写出高性能的 MATLAB 代码，主要参考官方文档中的提升性能的方法。

基础

关于代码结构

MATLAB 提供了很多方式来执行代码：命令行 vs 脚本 vs 函数，绝大多数代码在不同方式中执行都是等效的，但是考虑优化就不是一回事了。

在涉及 JIT 优化和作用域分析时，函数形式通常更有利于 MATLAB 对代码进行分析，尤其是作用域与类型推断，进而对代码运行进行优化，因此运行效率更高；命令行和脚本则更偏向交互式使用。

在命令行中交互式地调用某些内置函数，相比于从脚本文件中调用，可能会有不一样的结果。

基于模块化编程的思维，避免使用过大的单一文件，考虑将其中重复使用的功能拆分为简洁的函数文件，这种做法可以降低首次运行的成本。

优先使用局部函数而非嵌套函数，嵌套函数可以直接访问外层函数的变量（按照引用捕获），这可能限制 MATLAB 对其进行的优化，更好的做法是将所有需要的变量以函数参数的形式显式传递。

关于内置函数

避免重载内置函数，尤其不要对标准 MATLAB 数据类重载内置函数。

在可选的情况下，优先考虑使用 MATLAB 的内置函数来实现数组的常用运算，因为内置函数可能是直接使用底层语言实现的，至少也经过了大量的优化，比我们直接实现的效率肯定更高。

虽然 MATLAB 提供了很多基于字符串解析的代码执行/函数调用/变量修改功能，支持查询和修改当前工作区变量，但是它们会带来额外的开销，不利于代码优化，并且会降低代码的灵活性和可维护性。

例如：

  • 避免使用查询MATLAB工作区状态的函数，例如 which、whos、exist、dbstack 等具有显著动态性质的函数，运行时的自检会消耗大量计算资源。
  • 避免使用 eval、evalc、evalin 和 feval(fname) 等基于字符数组解析来执行的函数，从文本间接计算会显著降低运行效率，尽可能使用函数句柄代替。

关于变量

类型稳定：虽然MATLAB的变量没有固定的类型，但是仍然应该尽量避免在使用中改变类型和形状，应该直接创建新变量，而不是复用现有变量。

避免全局变量：全局变量会降低代码的可读性和性能，阻止代码优化，应该尽可能使用局部变量，通过显式传参消除全局变量。

避免数据硬编码：避免在代码中直接存储数据，例如硬编码大规模的字面量矩阵（例如超过 100 行的数据），将数据单独存储为 .mat 或 .csv 文件，在代码中通过命令读取文件来生成矩阵，无论从什么角度来说都是更好的做法。

好的代码主要从两个方面考虑：对人友好，无论是可读性还是可维护性等；对机器友好，便于运行和优化。上面提到的点既不是对人友好的，也不是对机器友好的。

数组优化

预分配+向量化

在使用数组时，需要考虑两个核心原则：预分配 + 向量化

  • 预分配数组：提前分配数组大小，避免在循环中动态扩展数组。（数组的动态扩容对性能的影响是很大的）
  • 向量化计算：尽可能使用矩阵和向量运算，避免 for 循环尤其是多层 for 循环。

我们可以对 for 循环中的几种常见写法进行简单的对比：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
clc;
clear;
close all;

num = 1e6;
freq = 0.01;

% case 1
clearvars -except num freq;
tic
x = sin(2*pi*(1:num)*freq);
toc

% case 2
clearvars -except num freq;
tic
x = zeros(1,num);
for k = 1 : num
   x(k) = sin(2*pi*k*freq);
end
toc

% case 3
clearvars -except num freq;
tic
for k = 1 : num
   x(k) = sin(2*pi*k*freq);
end
toc

% case 4
clearvars -except num freq;
tic
x = [];
for k = 1 : num / 5
   x = [x sin(2*pi*k*freq)];
end
toc

在我的电脑上测试结果如下

1
2
3
4
Elapsed time is 0.003326 seconds.
Elapsed time is 0.014658 seconds.
Elapsed time is 0.061899 seconds.
Elapsed time is 15.069429 seconds.

很显然，第一种完全向量化的写法是最优的，计算效率是最高的，
第二种写法对数组的整个空间进行了预分配，在for循环中对每一个元素赋值，效率虽然比不了向量化的写法，但是比后两种需要不断扩容数组的写法要快很多。

至于最后的两种写法：

  • 第三种写法性能比较低，因为每次都在向量中添加一个元素，尺寸在不断扩充，但是扩容都是就地进行的，由于预留空间的存在，并不是每次扩容都触发了整个数组的拷贝，因此实际效率并没有想象的那么糟糕。
  • 第四种写法的性能毫无疑问是最差的（这里已经将循环次数改小了），因为它是通过构造新矩阵并赋值的方式完成的矩阵扩容效果，在每一次 for 循环中，x = [x y] 都对整个数组进行了完整的拷贝，这导致计算复杂度更高，效率慢的可怕。

对 for 循环的厌恶和对向量化的推崇可能是 MATLAB 社区中的主流观点，但是 for 循环其实并没有那么糟糕，高版本的 MATLAB 可以对结构良好的 for 循环进行 JIT 优化。相比之下，数组的动态扩容和频繁拷贝才是真正糟糕的做法。

稀疏数组

稠密数组和稀疏数组在计算效率和内存开销上有巨大差异，如果一个矩阵确实是稀疏的，那么应该使用稀疏矩阵存储，
至少确保在主要计算环节中都将其存储为稀疏矩阵并处理，而不是转换为含有很多 0 的稠密矩阵。

例如对于稀疏矩阵来说，稀疏表示至少比稠密表示的内存占用少一个量级，如果还涉及到高维问题的离散，或者设计矩阵的 kron 运算，两者的存储和运算效率的差异更是非常巨大。

考虑下面的示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
clear; clc;

n = 15;
e = ones(n,1);
A = spdiags([-e 2*e -e], -1:1, n, n);

levels = 1:4;  % number of Kronecker levels
fprintf('Levels | Dense size (MB) | Sparse size (MB)\n');
fprintf('---------------------------------------------\n');

for k = levels
    Mf = 1;
    Ms = 1;
    for i = 1:k
        Mf = kron(Mf, full(A));
        Ms = kron(Ms, A);
    end

    % get memory info
    info_dense  = whos('Mf');
    info_sparse = whos('Ms');

    fprintf('%6d | %14.2f | %15.2f\n', k, info_dense.bytes/1024^2, info_sparse.bytes/1024^2);
end

这里从 15*15 的三对角矩阵开始，采用多层递归的 kron 运算构建大数组。运行结果如下

1
2
3
4
5
6
Levels | Dense size (MB) | Sparse size (MB)
---------------------------------------------
     1 |           0.00 |            0.00
     2 |           0.39 |            0.03
     3 |          86.90 |            1.24
     4 |       19553.30 |           52.55

可以看到，在稠密数组需要 19G 的存储时，稀疏版本仍然只需要 50 MB。

上述测试代码已经是 32GB 个人笔记本的极限：扩大 n 或者 k 都会导致内存不足的报错。

实践中使用大规模稀疏矩阵的情景通常对性能也非常敏感，稀疏矩阵的部分运算开销其实也有点大（例如频繁构造稀疏矩阵），也需要尽量避免：

  • 例如对于需要大量循环执行的代码，尽量把稀疏数组的构造放在循环之外。
  • 对于稀疏对角矩阵，可以考虑直接用向量存储，但是运算时要注意细节，因为 MATLAB 会自动对维度广播，可能写错。

尽可能使用稀疏矩阵和对应算法是一个非常普适的原则：

  • 例如在 C++ 中不要试图直接用嵌套的 std::vector<double> 存储大规模稀疏矩阵，必须手动实现稀疏矩阵的存储方式（专门为稀疏矩阵设计的数据结构），否则内存很可能撑不住，而且计算效率会非常低；
  • 不使用稀疏矩阵对应的数值算法还会直接影响数值结果的可靠性，尤其是大规模稀疏线性方程组的求解；
  • 对于 C++，建议的做法当然是使用现有的成熟第三方库，例如 Eigen。更好的做法是直接换一个语言，例如 MATLAB 或 Julia，C++ 从头开始写数值算法太难了。

并行计算

下面整理一些关于 MATLAB 并行计算的内容，主要是 Parallel Computing Toolbox 的内容，暂不涉及集群层面的并行，也不涉及 GPU 的使用。主要参考下列官方文档：

在涉及到并行计算时，存在很多新问题需要处理，例如其它进程并不会共享主进程的 persistent 变量等。
并行还会对通常的结果输出造成显著影响：（不是 BUG）

  • 本应向控制台输出的信息可能丢失，也可能乱序输出；
    • parfor 和 spmd 的标准输出可能会乱序
    • parfeval 的标准输出完全丢失
    • 建议向文件中输出，并且每个进程/线程打开不同的文件
  • 本应弹出的绘图窗口不再出现等。（建议直接保存图片并关闭 figure 对象）

parfor

将 for 循环改成 parfor 循环是最简单直接的并行方式，它会启动不同的进程/线程来同步执行循环，这对循环中的内容有一定的要求：

  • 每次迭代相互独立，结果不依赖于其他迭代
  • 不支持某些动态操作变量的函数，如 eval 和 assignin
  • 循环次数必须明确，不能依赖于运行时的计算结果
  • 不支持嵌套的 parfor 循环

对于简单的循环体，可以尝试直接用矢量化的语法，因为 parfor（以及其它的并行技术）在启动和结束过程中会带来额外的时间成本。

在 parfor 循环体内部，变量被大致分为如下几类：

  1. 循环变量，对它的限制是在循环内不能对循环变量再次赋值，并且不要在循环体外部再次使用，因为并行会导致无法确定循环变量的具体值
1
2
3
4
parfor i = 1:n
    i = i + 1; % not allowed
    a(i) = i;
end
  1. 切片变量，是指每个循环只访问该变量的特定位置，具体访问位置跟循环变量有关。在每一个循环中访问变量的位置必须是固定且不重合的。如果该变量在循环内被赋值，访问还必须是内存连续的（此时只能是循环变量再加固定的平移量），并且该变量不能在循环内动态变换大小
1
2
3
4
5
6
7
8
9
parfor i = 1:n
    x(i) = a(2*i*i);  % allowed;
    y(i+2) = a(i) + b(i+1); % allowed;

    c(i+1) = c(i) + 1; % not allowed;
    z(2*i) = i; % not allowed;
    a(i) = [];  % not allowed;
    a(end+1) = i; % not allowed
end
  1. 广播变量，是指纯粹的外部变量，在循环内未被重新赋值，只是作为只读变量被使用。
1
2
3
4
5
6
7
n = 100;
a = 10;
result = zeros(n,1);

parfor i = 1:n
    result(i) = i + a; % a 是广播变量
end
  1. 临时变量，指在循环内部创建的临时变量，该变量仅在单次循环中有效，在并行的不同循环中会创建不同的副本，不会相互影响，并且该变量不会通过索引变量引用（否则该被归类到切片变量）。
1
2
3
4
5
6
7
n = 100;
result = zeros(n,1);

parfor i = 1:n
    temp = i^2; % temp 是临时变量
    result(i) = temp;
end

除此之外，还有一类特殊的变量，它会在每次迭代中累加或减少，MATLAB 会自动处理这些变量，并在并行计算的循环结束后将修改汇总，以确保结果的正确性。

1
2
3
4
5
s = 0;

parfor i = 1:10
    s = s + i;
end

原则上来说，我们不应该在并行中执行这样的操作，在一般的并行框架下，这种数据汇总操作可能需要加锁才能保证安全性，但是 MATLAB 还是为其提供了有限支持。

parfeval

parfeval 允许我们在后台异步执行函数，并在计算完成时获取结果，这对于需要并行执行多个独立任务并在它们完成后收集结果的情况非常有用。（在后台执行意味着 worker 的标准输出会被丢弃）

对比 parfor 和 parfeval：

  • parfor 适合数据级的并行，parfeval 适合任务级的并行；
  • parfor 会阻塞执行流直到所有并行计算结束，parfeval 不会阻塞执行流，完全异步执行；

一个典型的例子如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
% 创建一个并行池
pool = gcp();

% 提交异步任务并获得 Future 对象
futures = cell(1, 10);
for i = 1:10
    futures{i} = parfeval(@myFunction, 1, i);
end

% 阻塞式等待所有任务完成并获取结果
results = zeros(1, 10);
for i = 1:10
    results(i) = fetchOutputs(futures{i});
end

% 显示结果
disp(results);

% 定义一个需要并行执行的函数
function result = myFunction(x)
    pause(2); % 模拟一个耗时操作
    result = x^2;
end

parfeval 函数的调用方式为 F = parfeval(fcn,numFcnOut,X1,...,Xm)，得到的 F 是一个 Future 对象，fcn 是需要调用的函数句柄，我们还需要提供返回值个数，以及所有的函数参数。这个语句会立刻返回，具体调用的函数会在后台继续异步执行。

我们可以通过 fetchOutputs 函数获取 Future 对象的返回值，这个语句会阻塞执行流，等待对应的后台任务完成并获取返回值。

我们可以使用 cancel 取消对应的任务

1
cancel(futures{i});

可以使用 wait 阻塞执行流，等待所有的后台任务完成

1
wait(futures);

parfevalOnAll 可以在所有的 worker 中异步执行相同的命令，通常是用于统一的环境配置或者清理工作，例如

1
parfevalOnAll(@clear,0,'mex');

这种语句如果直接放在普通的并行语句中可能会报错，就像在 parfor 中不允许调用 eval 一样，因此 MATLAB 专门提供了 parfevalOnAll。

与之类似的函数是 pctRunOnAll，与 parfevalOnAll 的最大不同是，它在 client 和 worker 之间是同步执行的，可能阻塞执行流。

注意：parfevalOnAll / pctRunOnAll 不仅在所有 worker 中执行，还会在 client 中执行。

此外还有如下的函数：

  • afterEach：在某一个 worker 完成时执行指定的回调函数
  • afterAll：在所有 worker 完成时执行指定的回调函数
  • fetchNext：获取当前已经完成的某个 worker 的结果

其中的 afterEach 需要注意，有多个同名但不同含义的接口（例如下面的 afterEach 函数是 parallel.pool.DataQueue.afterEach 的方法，控制消息队列收到信息的行为），只有在类型正确时才会调用正确的版本。

需要注意的是，因为 parfeval 会在后台执行任务，通常的控制台输出（包括警告信息）都会被隐藏。
可以使用 parallel.pool.DataQueue 收集各个 worker 的输出并在 client 端进行输出，例如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
% 创建一个并行池
pool = gcp();

% 在 client 创建一个 DataQueue
dq = parallel.pool.DataQueue;
% 添加回调函数，收到 msg 即打印
afterEach(dq, @(msg) fprintf('%s\n', msg));

% 提交异步任务并获得 Future 对象
futures = cell(1, 10);
for i = 1:10
    futures{i} = parfeval(@myFunction, 1, i, dq);
end

% 阻塞式等待所有任务完成并获取结果
results = zeros(1, 10);
for i = 1:10
    results(i) = fetchOutputs(futures{i});
end

% 显示结果
disp(results);

% 定义一个需要并行执行的函数，输出信息到 DataQueue
function result = myFunction(i, dq)
    % 发送消息到 client
    send(dq, sprintf('[task %d] started', i));

    pause(1); % 模拟一个耗时操作
    result = i^2;

    % 发送消息到 client
    send(dq, sprintf('[task %d] finished, result = %d', i, result));
end

spmd

spmd 是一种类似于 MPI 的并行实现方式，支持在不同 worker 之间进行数据交换：发送数据，接收数据，等待接收等，因此适合需要跨 worker 频繁通信的场景。

示例如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
spmd
    workerID = spmdIndex;

    % 每个worker执行其计算任务
    data = workerID^2;

    if workerID ~= 1
        % 其他的wokrer将其结果发送给worker1
        spmdSend(data, 1);
    else
        % wokrer1接收来自其他wokrer的结果并进行汇总
        total = data; % 自己的计算结果
        for i = 2:spmdSize
            % wokrer1接收其它worker的结果
            received_data = spmdReceive(i);
            fprintf('Worker 1 received data %d from worker %d\n', received_data, i);
            total = total + received_data; % 求和
        end
        fprintf('Total result is %d\n', total); % 输出汇总结果
    end
end

与 parfor 类似，spmd 语句块整体是阻塞式的，会阻塞执行流直到所有 worker 都完成计算。

与 parfor 和 parfeval 不同，spmd 不存在并行任务数目不等于 worker 数目的情况：始终都是在所有 worker 同时执行一份代码（只是 spmdIndex 等不同，用于相互通信），因此并行数目必然等于 worker 数目。

parpool

前面的所有并行语句其实都需要依赖并行池（parpool），MATLAB 在执行中遇到 parfor 等语句时，如果没有开启并行池，会自动尝试开启并行池，因此我们通常可以忽略并行池的存在，但是 MATLAB 实际上也提供了更精细地控制并行池的相关接口。

对于 R2025a 之前的版本: 一个 MATLAB 客户端会话中一次只能有一个并行池。

parpool 函数可以启动并行池，可以指定 worker 数量，也可以使用默认值

1
2
3
parpool();
parpool('local');
parpool('local', 4);

worker 数量的默认值通常是本地的核数（也可以通过配置修改），这也是允许的最大值，超过这个值会直接报错。

gcp函数可以获取当前并行池对象（并行池是一个全局的单例对象）

1
2
3
pool = gcp(); % 获取当前并行池对象，如果没有并行池，会自动启动并返回

pool = gcp('nocreate'); % 尝试访问但不创建并行池，如果没有并行池，则返回空

注意这里在无参数情况下调用 gcp() 时会自动启动并行池，如果仅仅需要查看状态必须加上 'nocreate' 参数。

我们可以通过下面的语句获取当前并行池的 worker 数量

1
2
3
4
5
6
poolobj = gcp('nocreate');
if isempty(poolobj)
    poolsize = 0;
else
    poolsize = poolobj.NumWorkers
end

下面的语句可以获取并手动关闭并行池

1
delete(gcp('nocreate'));

这是建议的写法，对空值进行 delete 没什么影响。

无参数形式则是错误的写法，因为如果并行池不存在，就会先启动再关闭。

1
delete(gcp());

我们也可以手动获取并保存 parpool 函数返回的并行池，然后通过 delete 函数删除。

1
2
3
4
5
poolobj = parpool();

...

delete(poolobj)

退出 MATLAB 时，未关闭的并行池也会随之自动关闭，因此无法用这种方式实现脱离 MATLAB 主进程后的任务长期运行。（这部分内容可以在 MATLAB 文档中检索 cluster，batch 等关键词，在使用中可能需要结合 PBS / Slurm / LSF 等任务调度系统）

MATLAB 的并行池支持多进程和多线程两种工作模式，通常使用的都是多进程模式，可以在设置界面中修改默认选项，在启动时也会展示并行池的基本信息

1
2
3
parpool('Processes');

parpool('Threads');

多线程是一个较新的工作模式，部分功能在多线程模式下受限或不支持，尤其是很多底层库（ BLAS / LAPACK / MKL）本身已经支持多线程，再叠加 thread pool 容易过度竞争。

多进程和多线程模式下的并行池区别大致就是多进程和多线程本身的区别。
可以使用下面的代码段检查当前并行池的类型

1
2
3
4
5
6
7
8
9
10
pool = gcp();

if isempty(pool)
    disp('No parallel pool is currently open.');
else
   disp(class(pool))
   % parallel.ThreadPool
   % or
   % parallel.ProcessPool
end

在 R2021b 版本之后，还有一个后台并行池 backgroundPool 的概念。

在实际计算中，我们需要关注并行池的 worker 数量，以及每一个 worker 可以使用的 maxNumCompThreads，这两个概念是非常重要的：

  • worker 太多会导致相互的计算资源冲突：每一个 worker 可以分配到的计算资源变少；
  • maxNumCompThreads 限制了在底层调用 BLAS 等函数自动进行多线程运算时可以使用的线程数，如果太大可能相互竞争，如果太小则会太慢。
    • MATLAB 的 client （主进程）默认使用的 maxNumCompThreads 等于物理核数；
    • MATLAB 在并行时， 每一个 worker 的 maxNumCompThreads 默认值为 1，这是非常关键的计算资源配置，对于 worker 数量很大的情况，maxNumCompThreads = 1 比较合理，但是在 worker 较少时，最好将其修改为 物理核数 / worker 数量。

我们还需要关注的是 parfor 和 parfeval 背后的调度规则：

  • 如果 worker 数量大于并行规模，多余的 worker 会被闲置；
  • 如果 worker 数量小于并行规模：
    • parfor 会根据 worker 数量和并行规模，自动将整个任务分块，将分块子任务分配给每一个 worker。这种处理对一些简单的并行计算是合适的，但是对于复杂的情况，这显然不能保证 worker 负载均衡，因此 MATLAB 还会进行有限的动态调度，例如给已经完成任务的 worker 分配剩余的任务。
    • parfeval 采用标准的任务队列模型，每一个 worker 获取一个任务执行，执行完成后再从剩余任务队列中获取下一个任务执行。

补充

哪些代码不应该追求性能

性能优化并不是无条件的编程目标，在很多实际场景中，过早或过度的性能优化不仅收益有限，反而会显著降低代码的可读性、可维护性和可复现性，还会阻碍代码调试。

下面是几类典型的，不需要追求性能的代码：

  • setup：包括一些必要的启动配置，参数解析
  • （静态）几何网格的读取/构造（如果使用是动网格/贴体网格，需要频繁地计算和更新网格，那就完全不一样了，这毫无疑问是性能瓶颈之一）
  • 数据加载与预处理
  • 输出目录、日志系统的配置
  • 计算结果输出

这些代码的执行次数很少，或者是受限于 I/O，并不需要刻意追求代码本身的高性能，而是应该追求代码的可读性和正确性。

一个合理的流程应该是：

  • 先写出清晰、正确的代码
  • 再通过 profiling 确认性能瓶颈
  • 仅对性能瓶颈部分进行针对性的优化

在进行 profiling 时，最好不要使用并行计算，并行可能会干扰性能测试的结果。

资源与性能检测

在实际使用时可能需要关注当前环境下的计算资源：

  • MATLAB 版本
  • 是否具有 GUI
  • 物理核心数（默认是并行池允许的最大数目）
  • client maxNumCompThreads
  • worker maxNumCompThreads
  • BLAS 版本
  • 系统内存信息

下面的测试脚本可用于检测这些关键信息

check.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
%% MATLAB check script

fprintf('=================================================\n');
fprintf('MATLAB version: %s\n', version('-release'));
fprintf('UseDesktop: %d\n', usejava('desktop'));
fprintf('Physical cores: %d\n', feature('numcores'));

clientThreads = maxNumCompThreads();
fprintf('Client maxNumCompThreads: %d\n', clientThreads);

% ---- Parallel pool ----
pool = gcp('nocreate');

if isempty(pool)
    fprintf('Parpool: none\n');
    fprintf('  Hint: run parpool("local") to enable parallel execution.\n');
    fprintf('  Hint: check client thread usage before benchmarking.\n');
else
    fprintf('Parpool: %s, %d workers\n', class(pool), pool.NumWorkers);

    spmd
        t = maxNumCompThreads;
    end

    workerThreads = unique([t{:}]);
    fprintf('Worker default maxNumCompThreads: %s\n', mat2str(workerThreads));

    if numel(workerThreads) > 1
        fprintf('  Warning: inconsistent maxNumCompThreads across workers.\n');
    end

end

fprintf('BLAS version: %s\n', strtrim(version('-blas')));

if ispc
    [~, mem] = system('systeminfo | findstr "Total Physical Memory"');
    fprintf('Memory (Windows):\n%s\n', strtrim(mem));

elseif isunix
    [status, cmdout] = system('free -h');
    fprintf('Memory (Linux):\n%s\n', strtrim(cmdout));

end

fprintf('=================================================\n');