如何对Matlab中的目标代码进行向量化以提升计算效率？

核心计算逻辑拆解 #

原代码中n_j的计算本质是对每个时间点b和索引j，计算以下标量：

n_j(b,j) = sum( exp(-1i * t(b) * delta_E) .* multi_i .* (V(j,:)' * V(j,:)), 'all' )

其中：

• delta_E(p,q) = E(p) - E(q)（L×L矩阵）
• multi_i(p,q) = alpha_i(p) * i_beta(q)（L×L矩阵，由两个向量外积得到）
• V(j,:)' * V(j,:)是矩阵V第j行的外积（L×L矩阵）

将求和式转化为矩阵运算形式，即可实现完全向量化，彻底消除循环。

向量化实现代码 #

直接替换原有的parfor和内层for循环，使用以下代码：

% 预计算时间相关的指数项：num×L×L的页数组，每个页对应一个t(b)的exp矩阵
exp_term = exp( -1i * reshape(t, num, 1, 1) * delta_E );

% 构造M数组：num×L×L，每个页M(b,:,:) = multi_i .* exp_term(b,:,:)
M = multi_i .* exp_term;

% 页矩阵乘法：计算V * M * V'，得到num×L×L的页数组
VMVt = pagemtimes( pagemtimes(V, M), V' );

% 提取每个页的对角线元素，得到num×L的n_j矩阵
n_j = diag(VMVt, 0, 3);

代码细节解释 #

1. 指数项的向量化构造

   • 将num×1的t重塑为num×1×1的三维数组，借助Matlab广播机制，与L×L的delta_E相乘后自动扩展为num×L×L的页数组exp_term。每个页对应一个时间点t(b)的exp(-1i*t(b)*delta_E)矩阵，彻底替代对b的循环。

2. 构造M数组

   • L×L的multi_i与num×L×L的exp_term逐元素相乘，广播后生成num×L×L的M数组，每个页对应单个时间点的权重矩阵。

3. 页矩阵乘法

   • 用pagemtimes处理页数组的矩阵乘法：先对每个页执行V * M(b,:,:)，得到num×L×L的中间结果；再与V'做页乘法，得到num×L×L的VMVt，每个页对应V*M(b,:,:)*V'矩阵。

4. 提取对角线元素

   • 通过diag(VMVt, 0, 3)在页维度上提取每个页的主对角线元素，直接得到num×L的n_j矩阵，与原代码循环计算的结果完全一致。

性能优势 #

• 完全消除双层循环，利用Matlab底层BLAS/LAPACK优化的矩阵运算，效率远高于循环实现；
• 避免parpool的启动开销和线程调度成本，在L和num较大时，性能提升尤为显著。

内容的提问来源于stack exchange，提问作者Mateusz Lisiecki