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

基础

关于代码结构

MATLAB 提供了很多方式来执行代码：命令行 vs 脚本 vs 函数，绝大多数代码在不同方式中执行都是等效的，但是考虑优化就不是一回事了，通常的运行效率关系为：命令行 < 脚本 < 函数，也就是说函数的执行速度通常是最快的。

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

优先使用局部函数而非嵌套函数，嵌套函数可以直接访问外层函数的变量（按照引用捕获），这会影响效率，更好的做法是将所有需要的变量以函数参数的形式显式传递。

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

关于内置函数

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

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

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

例如：

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

关于变量

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

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

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

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

数组优化

预分配+向量化

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

  • 预分配数组：提前分配数组大小，避免在循环中动态扩展数组。
  • 向量化计算：尽可能使用矩阵和向量运算，避免 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] 都对整个数组进行了完整的拷贝，这导致计算复杂度更高，效率慢的可怕。

稀疏数组

稠密数组和稀疏数组在计算效率和内存开销上有巨大差异，如果一个矩阵确实是稀疏的，那么应该使用稀疏矩阵存储， 至少确保在计算的主要环节中都将其存储为稀疏矩阵并处理，而不是转换为含有很多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 存储稀疏矩阵，必须手动实现稀疏矩阵的存储方式，否则内存很容易撑不住，而且计算效率也会非常低。当然，更好的做法是用现有的成熟第三方库，例如 Eigen 的稀疏数组。

并行计算

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

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允许我们在后台异步执行函数，并在计算完成时获取结果，这对于需要并行执行多个独立任务并在它们完成后收集结果的情况非常有用。

一个典型的例子如下 

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。

此外还有如下的函数：

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

其中的afterEach很奇怪，有两个同名但不同含义的接口，只有在类型正确时才会调用正确的版本。 这些函数的使用情景较少，并且并不是必要的，因此省略。

需要注意的是，因为parfeval会在后台执行任务，通常的控制台输出（包括警告信息）都会被隐藏。

spmd

与parfor和不同，spmd是一种类似于MPI的并行实现方式，除了支持相互独立的循环并行，还支持在不同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

parpool

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

parpool函数可以启动并行池 

1
2
parpool(); % 启动本地默认配置的并行池
parpool('local', 4); % 启动具有4个工作进程的本地并行池

这里无参数情况下启动的本地默认配置与本地的核数有关，通常是尽可能开到最大。

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

1
2
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')); % 关闭当前并行池

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

1
2
3
4
5
poolobj = parpool();

...

delete(poolobj)

在MATLAB退出之后，未关闭的并行池也会随之自动关闭。

MATLAB的并行池其实支持多进程和多线程两种工作模式，可以在设置界面中修改默认选项，在启动时也会展示相应的信息 

1
2
3
parpool("Processes"); % 启动多进程并行池

parpool("Threads"); % 启动多线程并行池

多进程和多线程的并行池在 MATLAB 中的对比如下表

特性 多进程并行池 (Processes) 多线程并行池 (Threads)
内存空间 独立内存空间，每个进程有自己的内存 共享内存空间，所有线程共享内存
启动和管理开销 高 低
通信开销 高 低
数据竞争 无数据竞争风险 存在数据竞争风险
资源消耗 较高，因独立内存占用较多 较低，共享内存占用较少
适用任务 计算密集型、需要隔离的任务 需要频繁数据传输、共享数据的任务
默认配置 是（默认是多进程） 否
隔离性 高 低

可以使用下面的代码段检查当前并行池的类型 

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

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

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

补充

资源与性能检测

在实际使用时需要关注MATLAB进行大规模科学计算所需要的计算资源：

  • MATLAB版本
  • 系统内存大小
  • 系统核心数（支持的并行数目）

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

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
%% Check MATLAB Environment

% Get current MATLAB version
matlabVersion = version('-release');
disp(['MATLAB Version: ', matlabVersion]);

% memory
if ispc
    % Windows
    [~, cmdout] = system('systeminfo | findstr "Total Physical Memory"');
    disp('Memory information (Windows):');
    disp(cmdout);
elseif isunix
    % Unix-like
    [status, cmdout] = system('free -h');
    disp('Memory information (Linux):');
    disp(cmdout);
else
    disp('System memory information not available on this platform.');
end

% Get maximum supported workers
maxWorkers = feature('numcores'); % Get maximum number of cores

% Check current parallel pool status, open if not already open, and execute test task
pool = gcp('nocreate'); % If no pool exists, do not create new
if isempty(pool)
    disp(['No parallel pool found. Opening a new one with ', num2str(maxWorkers), ' workers...']);
    parpool('local', maxWorkers); % Open local pool based on maximum supported cores
    disp('Parallel pool opened successfully.');
else
    disp(['Current parallel pool size: ', num2str(pool.NumWorkers)]);
end

% Execute parallel task to output "Hello, world!" using maximum supported workers
parfor idx = 1:maxWorkers
    fprintf('Hello, world! (from %d)\n', idx);
end

%  close parallel pool
delete(gcp);