提高matlab运算效率——预分配内存

提高 matlab 的运算效率有很多技巧，比如用向量运算代替循环运算等。
这次介绍最常见的提高运算效率的方法——预分配内存。
预分配内存是一件看上去非常平常以至于经常无视的技巧，但面对大量循环运算时，预分配内存能节约的时间是非常可观的。

最浪费时间的写法

下面这种最浪费时间的写法竟然是自己最常使用的，可见平时编程习惯很欠优化。

clear
tic
a_array = [];
for i = 1:1e5
    a = sin(i);
    a_array = [a_array;a];
end
toc

上面这种不断在原数组后补充新元素组成新数组的做法十分浪费时间。

• 循环 1e4 次耗时 0.020635 秒。
• 循环 1e5 次耗时 1.427629 秒。是 1e4 次的 69.18 倍。 1e4 次的 69.18 倍。
• 循环 1e6 次耗时 763.709362 秒。是 1e5 次的 534.95 倍！

很明显看到，耗时并不随着循环次数线性增长，即并不是循环次数变为 10 倍则运行耗时也变为 10 倍。可见循环次数越多，这种编程方式的缺陷就越严重。

另一种没有预分配内存的写法

下面这种写法是按索引不断增补元素，运行效率比前一种有明显提高。但这种写法和预分配内存相比还是效率太低了。

tic
for i = 1:1e6
    b = sin(i);
    b_array(i) = b;
end
toc

tic
for i = 1:1e6
    b = sin(i);
    b_array(i) = b;
end
toc

上面的程序运行结果为

• 历时 0.092180 秒。
• 历时 0.018310 秒。
  可以看到相同的程序第二次运行耗时比一次要短得多。这是因为程序在第一次运行时已经为 b_array 分配了内存，因此节省了时间。由于不需要分配内存的时间，甚至可以看到这会比下面预分配内存的写法运行耗时还要短，这也侧面印证了预分配内存的有效性。

从上面的统计结果来看运行效率远远高于第一种方法。

预分配内存的写法

clear
tic
c_array = zeros(1, 1e5);
for i = 1:length(c_array)
    c = sin(i);
    c_array(i) = c;
end
toc

因此在实际编程中应为数组预分配内存，如果数组长度预先无法知晓，可以先预估一个长度分配。但是需要注意，如果预估的偏长则需要在最后将数组截断，否则数组后面会挂着一串零。
如果实在没有办法预估，也可以用第二种按索引增补元素的方法。但是如果增补的是数组时，那么还是采用预分配内存的方式更明智，具体原因见下。

面对补充数组时

当每次补充的是数组而不是单一元素时，预分配内存的优势将更加明显。
其实这种场景要更加常见，因为动态系统仿真时状态量往往不止一维，即状态量往往是数组。
此时对比上面第 2 和第 3 种写法的运算效率如下

%% 按索引补充数组
clear
tic
for i = 1:1e6
    b = [sin(i),sin(2*i),sin(3*i),sin(4*i),sin(5*i),sin(6*i),sin(7*i)];
    b_array(i,:) = b;
end
toc

%% 预分配内存
clear
tic
c_array = zeros(1e6,7);
for i = 1:length(c_array)
    c = [sin(i),sin(2*i),sin(3*i),sin(4*i),sin(5*i),sin(6*i),sin(7*i)];
    c_array(i,:) = c;
end
toc

可以看出按索引增补数组的方法运算耗时呈爆炸式增长，在循环运行 1e6 次时耗时已经不可接受了。

面对补充数组时还有另一个常见的情况即数组的转置。这在系统的动态仿真中也是十分常见的。因为系统状态量一般习惯写为列向量，而最终结果为了方便和时间列向量对应又希望状态矩阵中每个时刻的状态是行向量，即每一行对应不同时刻。这就难免涉及到数组的转置。
对比下面几种写法的运行效率

clear
tic
c_array = zeros(1e8,7);
for i = 1:length(c_array)
    c = [sin(i),sin(2*i),sin(3*i),sin(4*i),sin(5*i),sin(6*i),sin(7*i)];
    c_array(i,:) = c;
end
toc

clear
tic
c_array = zeros(1e8,7);
for i = 1:length(c_array)
    c = [sin(i);sin(2*i);sin(3*i);sin(4*i);sin(5*i);sin(6*i);sin(7*i)];
    c_array(i,:) = c';
end
toc

clear
tic
c_array = zeros(7,1e8);
for i = 1:length(c_array)
    c = [sin(i);sin(2*i);sin(3*i);sin(4*i);sin(5*i);sin(6*i);sin(7*i)];
    c_array(:,i) = c;
end
c_array = c_array';
toc

最终运行结果为

• 第 1 种运行耗时 113.274174 秒
• 第 2 种运行耗时 157.922921 秒
• 第 3 种运行耗时 129.031270 秒
  其实这里第 1 种和第 2 种在实际使用时是差不多的。要么就是在一开始就把状态向量都写成行向量，但这样在各个函数内部会很不舒服；要么就是在保存时把列向量转置成行向量。
  是实际使用时采用第 2 种还是第 3 种写法都可以，差别没有那么大。