如何加速MATLAB代码运行

学习笔记

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=4096;
theta=[1    1    2    3    1    1    2    3    1    2    3    0    2    3    0    2    3    0    1    3    0    1    3    0    1    2];
fai=[1787    1077    1753    697    1523    5    2035    331    1920    130    4    85    551    15    1780    1960    3    145    1019    691    132    42    393    502    201    1064
1502    602    749    1662    1371    9    131    1884    1268    1784    19    1839    81    2031    76    336    529    74    68    186    905    1751    1516    1285    1597    1712
1887    521    590    1775    1738    2032    2047    85    1572    78    26    298    1177    1950    1806    128    1855    129    269    1614    1467    1533    925    1886    2046    1167
1291    301    1353    1405    997    2032    11    1995    623    73    1839    2003    2019    1841    167    1087    2032    388    1385    885    707    1272    7    1534    1965    588];
A = zeros(M);
B = eye(M);
L = 0:M-1;
for matrixNum = 1:14
    t_k = theta(matrixNum);
    f_4i_M = floor(4*L/M);
    f_k = fai(f_4i_M+1,matrixNum)‘;
    col_1 = M/4*(mod((t_k+f_4i_M),4)) + ...
        mod((f_k+L),M/4);
    row_col = col_1+1 + L*M;
    C_temp = zeros(M);
    C_temp(ind2sub([M,M],row_col)) = 1;
    C{matrixNum} = sparse(C_temp)‘;
end
H = [A A A B B+C{1};B+C{8} B+C{7}+C{6} A A B;A B B+C{5} A C{4}+C{3}+C{2}];
H_23 = [A A;B C{11}+C{10}+C{9};C{14}+C{13}+C{12} B];
H=[H_23 H];
H_full = full(H);
whos H H_full
%% 稀疏矩阵乘法测试
message = randi([0 1],1,28672);
tic;H*message‘;toc;
tic;H_full*message‘;toc;

%% 数据类型测试
Gc = randn(16384);
Gc_logic = Gc>0;
a=randi([0 1],1,16384);
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([0 1],1,1000000);

tic;
encodeData_extend = [1 encodeData];
for num = 2:length(encodeData_extend)
    encodeData_extend(num) = xor(encodeData_extend(num),encodeData_extend(num-1));
end
toc;
tic;
encodeData1 = [1 encodeData];
encodeData1_sum = cumsum(encodeData1);
encodeData_2 = mod(encodeData1_sum,2);
toc;