如何加速MATLAB代码运行

2023-01-07 19:54:27

学习笔记

V1.0 2015/4/17

如何加速MATLAB代码运行

概述	本文源于LDPCC的MATLAB代码，即《CCSDS标准的LDPC编译码仿真》。由于代码的问题，在信息位长度很长（大于10000）情况下，代码无法正常运行或执行速度很慢。本文将叙述代码修改过程中的一系列手段，然对其加速原理不做探究

修订历史

以下表格展示了本文档的修订过程

日期	版本号	修订内容
2015/04/17	V1.0	初始版本

简介

本程序基于MATLAB 2014a 编写，本文档中提到的"MATLAB"均指该特定版本MATLAB。代码运行结果测试机器是T530i i3 3110M 16G @1600MHz。

MATLAB的帮助文档 - Advancde Software Development - Performance and Memory中提及了一些改善代码性能的一些手段。比较通用的包括向量的预先分配内存，这一点在编辑器里也会提示。有时候预先分配内存与否和性能关系很大，譬如

tic

x = 0;

for i = 2:1000000

x(i) = x(i-1)+5;

end

toc

tic

x = zeros(1,1000000);

for i = 2:1000000

x(i) = x(i-1)+5;

end

toc

运行结果显示为"时间已过 0.500985 秒"和"时间已过 0.073622 秒"。另外在声明变量的时候不使用原有变量，而创建新变量也可以减少运行时间。还有包括选择'&'和'&&'的差异。

MATLAB还提供了一些改善性能的手段，包括

将长脚本拆开成小段，调用执行；
将大的代码块分开为独立的函数；
将过分复杂的函数或是表达式采用简单的来代替；
采用函数，而不是脚本；
向量化代码，采用MATLAB自带的函数；
采用矩阵的稀疏结构；
运行MATLAB的时候不要在后台运行其他大的程序；
不要重载任何MATLAB的内建函数或数据类型。

不得不说上面的技巧有很多是废话，而其中向量化是最有效的一种方法之一。向量化代码中有很多常用的函数，包括

函数	功能（暂略）
all
any
cumsum
diff
find
ind2sub
ipermute
logical
meshgrid
ndgrid
permute
prod
repmat
reshape
shiftdim
sort
squeeze
sub2ind
sum

在代码撰写修改过程中，可以多考虑考虑以上函数。

实例

实例来自于《CCSDS标准的LDPC编译码仿真》中代码（实际上有点点差别），代码优化从以下几个方面进行

稀疏
类型转换
向量化

稀疏

仿真中的第一个困难在于ccsdscheckmatrix函数在输入SIZE_M很大的时候，先不说运行时间，直接就爆内存了。（输入参数4096，2/3）

先分析分析内存的问题，实际上这个函数的最后输出结果就是一个矩阵，这个矩阵的大小是12288×28672，计算double型的内存占用也就2G左右。但是函数运行过程中产生了很多中间变量没有清除。当然最后的解决办法也没有去管这些东西，由于矩阵H是稀疏矩阵，所以之际采用sparse后，这个运行就没有任何问题了。

对于矩阵H和H_sparse = spares(H)，占用内存如下（当然H要是稀疏的，不然得不偿失）

Name 　　　　Size 　　　　　　　　Bytes 　　　　　Class　　 Attributes

H 　　　　　　12288x28672 　　2818572288 　　double

H_sparse　 12288x28672　　 1736712 　　　　double 　　sparse

也可以对比稀疏矩阵和原始矩阵的运行时间（和稀疏程度有关）

代码：tic;H*message';toc;

结果：时间已过 0.288934 秒。

代码：tic;H_sparse*message';toc;

结果：时间已过 0.001210 秒。

类型转换

MATLAB中的运算符支持多种类型，譬如矩阵乘法中多用double型变量，但如果一个矩阵是逻辑输入也没有关系。但运算速度差异较大，譬如

>> Gc_logic = Gc>0;

>> a=randi([0 1],1,16384);

>> tic;b = a*Gc;toc

时间已过 0.107618 秒。

>> tic;b = a*Gc_logic;toc

时间已过 0.503132 秒。

观测结果类型为double，我们可以大胆推测实际上逻辑型变量在运算过程中先转化为了double型（逻辑怎么乘呢？）另一个实验结果是

>> tic;Gc_logic=double(Gc_logic);b = a*Gc_logic;toc

时间已过 0.546412 秒。

这一定程度上证明了我们的假设。所以在运算过程中数据类型是重要的，如果上述乘法出现在循环内，那么实现转化矩阵类型是必要的。即使只运行一次，那么显式的转化矩阵类型（特制新建变量）也有好处。譬如

>> tic;Gt=double(Gc_logic);b = a*Gt;toc

时间已过 0.373506 秒。

通过创建新变量，运行速度些许。

向量化

向量化实际上是原代码修改中获益最大的方法，这实际上是因为原先的译码程序写了太多的循环。向量化后运行时间变成了原先的1/40 。当然，原先的代码通用性强，而向量化这个过程实际上是运用了H的一些结构的。译码函数太复杂，此处不做举例。

此处分析差分调制中的例子（实际上对这个程序没有什么影响）

原来的代码是这个样子的（更新值为其本身和前一个值的异或）

encodeData_extend = [1 encodeData];

for num = 2:length(encodeData_extend)

encodeData_extend(num) = xor(encodeData_extend(num),encodeData_extend(num-1));

end

向量化的结果为（累加模二代替异或）

encodeData1 = [1 encodeData];

encodeData1_sum = cumsum(encodeData1);

encodeData_2 = mod(encodeData1_sum,2);

运行时间分别为

时间已过 0.023424 秒。

时间已过 0.015003 秒。

虽然后者没有快很多，但这取决于向量的长度，长度大的话会有较大差距。

其他

MATLAB中提及的都能对代码运行速度带来细微的改进，包括

将长脚本拆开成小段，调用执行；
将大的代码块分开为独立的函数；
将过分复杂的函数或是表达式采用简单的来代替；
采用函数，而不是脚本；

上述测试脚本（和以上运行条件有差别）

%% 稀疏矩阵测试

M=;

theta=[                                                                                                    ];

fai=[                                                                                                    

                                                                                                    ];

A = zeros(M);

B = eye(M);

L = :M-;

for matrixNum = :

    t_k = theta(matrixNum);

    f_4i_M = floor(*L/M);

    f_k = fai(f_4i_M+,matrixNum)';

    col_1 = M/*(mod((t_k+f_4i_M),)) + ...

        mod((f_k+L),M/);

    row_col = col_1+ + L*M;

    C_temp = zeros(M);

    C_temp(ind2sub([M,M],row_col)) = ;

    C{matrixNum} = sparse(C_temp)';

end

H = [A A A B B+C{};B+C{} B+C{}+C{} A A B;A B B+C{} A C{}+C{}+C{}];

H_23 = [A A;B C{}+C{}+C{};C{}+C{}+C{} B];

H=[H_23 H];

H_full = full(H);

whos H H_full

%% 稀疏矩阵乘法测试

message = randi([ ],,);

tic;H*message';toc;

tic;H_full*message';toc;

%% 数据类型测试

Gc = randn();

Gc_logic = Gc>;

a=randi([ ],,);

tic;b = a*Gc;toc

tic;b = a*Gc_logic;toc %逻辑型运行花费时间

tic;Gc_logic=double(Gc_logic);b = a*Gc_logic;toc %类型转换

tic;Gt=double(Gc_logic);b = a*Gt;toc %建立新变量

%% 向量化测试

encodeData = randi([ ],,);

tic;

encodeData_extend = [ encodeData];

for num = :length(encodeData_extend)

    encodeData_extend(num) = xor(encodeData_extend(num),encodeData_extend(num-));

end

toc;

tic;

encodeData1 = [ encodeData];

encodeData1_sum = cumsum(encodeData1);

encodeData_2 = mod(encodeData1_sum,);

toc;