Matlab高效编程：向量化(vectorization)、矩阵化、变量预定义

笨牛慢耕

已于 2022-11-02 08:36:54 修改

阅读量6.5k

点赞数 10

分类专栏： Matlab 文章标签： matlab 矩阵向量化 vectorization 傅里叶变换

于 2022-11-01 22:35:17 首次发布

本文链接：https://blog.csdn.net/chenxy_bwave/article/details/127624235

版权

Matlab 专栏收录该内容

23 篇文章 38 订阅

订阅专栏

0. 前言

本文介绍几个Matlab常用的提高运行效率的编程技巧。

对一个基于数值化的方式计算一个连续函数的频谱（傅里叶变换）的例子给出了三种实现（双重循环、向量化实现、矩阵化实现）代码，对比了运行时间。

1. 变量预定义

由于Matlab是解释性执行语言，Matlab不要求变量使用之前进行预定义，也不要求内存的预分配，一切都可以在执行过程中动态分配。这种灵活的方式方便了用户编程，但是同时也带来了潜在的低效率问题。

以下针对一个简单的计算例子，用几种不同的代码实现方式来比较不同的写法对于运行效率的影响。

代码中用tic 和 toc这两个matlab的基本的定时器函数，来进行运行时间统计。为了避免已分配内存对运行时间的影响，每次新的实现方式的运行之前都用“clearvars”进行内存清除，当然num和freq不需要清除，所以用-except选项将它们保留。

clearvars;
close all;
clc;

num = 1e6;
freq = 0.01; % Normalized frequency

% case1
tic
for k = 1 : num
    sin(2*pi*k*freq);
end
toc

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

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

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

运行结果如下（Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz，8GB内存，64比特操作系统）：

历时 0.204018 秒。
历时 0.223268 秒。
历时 0.076024 秒。
历时 2.713520 秒。

case1和case2（差异在于是否用运算结果对x(k)进行显式赋值）运行时间差异不大，这个虽然看上去略微有点奇怪，但是这只是说明了即便是只计算sin(2*pi*k*freq)不显式地对x(k)进行赋值也是需要动态地分配临时内存以存储计算结果的。

case3相比case1,2快了30倍左右，仅仅是预先对x进行了定义（相应地进行内存分配）而已，就能够有如此巨大的速度提升，可见变量预定义（尤其是大的数组变量）对于高效的Matlab编程有多么重要！

case4是为了对比给出的一个更为恶劣的例子（惭愧的是，笔者曾很长时间用这种方式实现动态数组增长，模拟python中的list.append()的行为。。。故此特地将它作为反面教材举出）。主要它在运行循环数是其它几种情况的十分之一的条件下，甚至比case1和case2都还要慢十倍以上，也就是说实际上慢100倍以上！

Last but not the least, 以上所说的变量预定义所带来的好处当然主要是指向量或者矩阵类型的变量。对于标量（scalar）类型变量则是否进行预定义没有什么影响。

2. 向量化，vectorization

向量化是另一个提高Matlab代码运行效率的主要手段。

与向量化处理相对应的是如上一节所示的for-loop的实现方式。以下用向量化实现方式实现以上相同的处理，直观地感受一下向量化实现对于matlab程序运行效率带来的提升效果。

clearvars;
close all;
clc;

num = 1e6;
freq = 0.01; % Normalized frequency

tic
k = 1 : num;
x = sin(2*pi*k*freq);
toc

运行结果：

历时 0.015811 秒。

运行时间仅为上一节最快的那种的五分之一，5倍速提升，值得你拥有。

3. 矩阵化

矩阵化其实是向量化从1维向2维的推广而已，本质上是一样的，都是并行处理优化。

矩阵化实现优化要求编程者将算法实现方式写成基于矩阵运算的方式，算法设计难度也提高了一些，对编程者的算法设计能力要求更高了。但是考虑到算法运行效率提高所带来的好处，这种算法设计技巧值得磨炼。

以下通过一个用matlab以离散方式近似计算连续信号傅里叶变换的实现来对比三种实现的运行效率（例子取自参考文献[1]）：(1) 用二重循环的方式；(2) 向量化 + for-loop；(3) 矩阵化。

例：求以下矩形脉冲的频谱（即傅里叶变换） $F(\omega), \ \omega \in [-8 \pi, 8 \pi]$ ：

$f(t) = \left\{\begin{matrix} 1 & |t|<\frac{1}{2}\\ 0 & otherwise \end{matrix}\right.$ (1)

根据定义有：

$F(\omega) = \int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt = \int_{t_1}^{t_2}f(t)e^{-j\omega t}dt$ (2)

注意，由于采用数值分析的方式进行近似计算，所以必须把积分范围限定在一个有限范围以内。以上假定 $[t_1, t_2]$ 代表函数f(t)的主要取值范围，并在这个区间内进行近似计算。

3.1 离散化

基于数值分析计数计算连续函数积分的第一步就是先将待求式进行离散化。

令 $T = t_2 - t_1$ ，假定将该区间划分为等间隔的N个区间，每个区间长度为 $\Delta t = \frac{T}{N}, \ t_2 = t_1 + N\cdot \Delta t$ .，进行待求积分离散化处理得到以下求和式（N越大，该求和式计算得到的结果就越接近待求式真正的结果）：

$F(\omega) = \sum\limits_{n=0}\limits^{N-1} f(t_1 + n\cdot \Delta t) e^{j \omega (t_1 + n\cdot \Delta t)} \frac{T}{N}$ (3)

进一步，还需要对频域进行离散化处理。假设 $[\omega_1, \omega_2]$ 为 $F(\omega)$ 的主要取值区间，并取其中K个离散点进行对应的采样值的计算，令 $\Delta \omega = \frac{\omega_2 - \omega_1}{K}$ ，则有：

$F(\omega_1 + k \Delta \omega) = \frac{T}{N} \sum\limits_{n=0}\limits^{N-1}f(t_1+n\Delta t)e^{-j (\omega_1 + k \Delta \omega) (t_1+n\Delta t)}$ (4)

3.2 双重循环实现

双重循环实现中，用内层循环实现对式(4)右边的累加，用外层循环实现对频域采样点的扫描，示例代码如下：

clearvars;
close all;
clc;

N = 2000;                         % number of time-domain sampling points
T = 20;                           % 时域区间长度
t = linspace(-T/2, T/2-T/N, N);   % 时域采样点
f = zeros(1,length(t));
f((t > -1/2) & (t < 1/2)) = 1;    % 定义时域信号

Omega = 16 * pi;                  % 频域区间长度
K     = 1000;                      % 频域采样点数
omg   = linspace(-Omega/2, Omega/2 - Omega/K, K); % 频域采样点
F1     = zeros(1,length(omg));

% 双重循环实现
tic
for k = 1:K
    for n = 1:N
        F1(k) = F1(k) + f(n) * exp(-j * omg(k) * t(n));
    end
    F1(k) = F1(k) * (T/N);
end
toc

3.3 向量化实现

以上式(4)显然可以看作是两个向量的内积（点积，inner product, dot product），将式(4)改写如下：

$F(\omega_1 + k \Delta \omega) = \frac{T}{N} \begin{bmatrix} f(t_1) & f(t_1+\Delta t) & \cdots & f(t_1+(N-1)\Delta t) \end{bmatrix} \begin{bmatrix} e^{-j (\omega_1 + k \Delta \omega) (t_1)}\\ e^{-j (\omega_1 + k \Delta \omega) (t_1+\Delta t)}\\ \cdots \\ e^{-j (\omega_1 + k \Delta \omega) (t_1+(N-1)\Delta t)} \end{bmatrix}$

(5)

基于这个表达式可以得到向量化实现，代码如下所示：

% 向量化实现
tic
for k = 1:K
    F2(k) = (T/N) * ( f * exp(-j * omg(k) * t).' );
end
toc

3.4 矩阵化实现

对式(5)可以进一步进行矩阵化处理可以得到：

$\begin{bmatrix} F(\omega_1 + 0 \Delta \omega)\\ F(\omega_1 + 1 \Delta \omega)\\ \cdots\\ F(\omega_1 + (K-1) \Delta \omega) \end{bmatrix} \\= \frac{T}{N} \begin{bmatrix} f(t_1) & f(t_1+\Delta t) & \cdots & f(t_1+(N-1)\Delta t) \end{bmatrix} \cdot \\ \begin{bmatrix} e^{-j \omega_1 t_1 } &e^{-j (\omega_1 + \Delta \omega) (t_1)} &\cdots &e^{-j (\omega_1 + (K-1) \Delta \omega) (t_1)} \\ e^{-j \omega_1 (t_1+\Delta t)} & e^{-j (\omega_1 + \Delta \omega) (t_1+\Delta t)} & \cdots & e^{-j (\omega_1 + (K-1) \Delta \omega) (t_1+\Delta t)}\\ \cdots & \cdots & \cdots & \\ e^{-j \omega_1 (t_1+(N-1)\Delta t) } & e^{-j (\omega_1 + \Delta \omega) (t_1+(N-1)\Delta t)} &\cdots & e^{-j (\omega_1 + (K-1) \Delta \omega) (t_1+(N-1)\Delta t)} \end{bmatrix}$ （6）

代码实现如下：

% 矩阵化实现
tic
exp_mat = exp(-j*(kron(t',omg)));
F3 = (T/N) * ( f * exp_mat );
toc

以上用kron()先计算出式(6)最右边的矩阵。关于张量积（克罗内克积、直积），有兴趣的小伙伴可以参考量子笔记：张量、张量的阶数与向量的维数、外积、张量积，此处不再赘述。

运行结果：

历时 1.192834 秒。
历时 0.071625 秒。
历时 0.062656 秒。

可以看到向量化实现相比双重循环实现有17倍的速度提升，相当可观。但是矩阵化相比向量化的提升的效果非常有限，但是矩阵化的代码最为简洁，也有可取之处。

可以进一步用以下assert语句检查以上三个结果是否完全一致：

assert(isempty(find(abs(F1-F2) > 1e-6)))
assert(isempty(find(abs(F1-F3) > 1e-6)))

当然也可以用“max(abs(F1-F2))”查询它们相互之间差异最大的元素的绝对值是多少。我的运行结果式10的负十六次方左右，大家可以自行检验。

【参考文献】

[1] 谷源涛等，信号与系统--Matlab综合实验，高等教育出版社

笨牛慢耕

关注

10
点赞
踩
36

收藏

觉得还不错? 一键收藏
打赏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录