如何高效实现矩阵按列对应相乘并生成结果向量？

现有两个尺寸相同的矩阵A和B，二者均为独立列向量的集合。需要对A与B的对应列执行矩阵乘法（*或mtimes），生成元素数等于列数的向量。循环实现虽简单，但因需执行数千次且矩阵尺寸≥2000×16~200，需采用向量化提升效率。此外还涉及权重矩阵W。

示例代码如下：

nDof  = 250;    % 行数
nCols = 16;     % 列数

% 示例复矩阵
A = rand(nDof,nCols) + 1j*rand(nDof,nCols); 
B = rand(nDof,nCols) + 1j*rand(nDof,nCols); 
% 权重矩阵
W = rand(nDof,nDof)  + 1j*rand(nDof,nDof); 

% 方案1：循环实现 - 逐列执行带/不带权重的矩阵乘法
X1 = zeros(nCols,1);  Y1 = zeros(nCols,1);
for ii = 1:nCols
    X1(ii) = A(:,ii)'     * B(:,ii);
    Y1(ii) = A(:,ii)' * W * B(:,ii);
end

% 方案2：尝试避免循环
X2temp = A'     * B;
Y2temp = A' * W * B;
% 取对角元素作为结果
X2 = diag(X2temp);
Y2 = diag(Y2temp);
    
% 验证结果差异
disp(X2-X1)
disp(Y2-Y1)

尝试通过直接矩阵乘法避免循环后，结果矩阵的对角元素与循环结果接近但不相等，误差超出浮点精度范围。测试结果如下：

X2-X1  =  1.0e-13 *

  0.426325641456060 + 0.355271367880050i
  0.284217094304040 + 0.852651282912120i
 -0.284217094304040 - 0.426325641456060i
  0.000000000000000 + 0.142108547152020i
 -0.284217094304040 - 0.213162820728030i
  0.284217094304040 - 0.923705556488130i
  0.568434188608080 + 0.355271367880050i
  0.000000000000000 + 0.142108547152020i
  0.000000000000000 + 0.071054273576010i
 -0.284217094304040 + 0.497379915032070i
  0.284217094304040 - 0.426325641456060i
 -0.284217094304040 + 0.142108547152020i
  0.000000000000000 + 0.000000000000000i
 -0.284217094304040 + 0.142108547152020i
 -0.284217094304040 - 0.213162820728030i
  0.284217094304040 - 0.284217094304040i

Y2-Y1 = 1.0e-10 *

  0.000000000000000 + 0.036379788070917i
 -0.018189894035459 + 0.072759576141834i
 -0.072759576141834 - 0.072759576141834i
  0.036379788070917 - 0.109139364212751i
  0.090949470177293 + 0.054569682106376i
 -0.109139364212751 + 0.054569682106376i
  0.109139364212751 - 0.163709046319127i
  0.018189894035459 - 0.127329258248210i
  0.072759576141834 - 0.018189894035459i
  0.036379788070917 - 0.036379788070917i
  0.018189894035459 + 0.000000000000000i
  0.109139364212751 - 0.072759576141834i
 -0.018189894035459 - 0.018189894035459i
  0.000000000000000 - 0.054569682106376i
 -0.018189894035459 + 0.036379788070917i
  0.036379788070917 - 0.036379788070917i

现寻求：正确的实现方式是什么？最快的实现方式是什么？

误差原因 #

你用A'*B取对角的方式，本质是计算了所有列对的内积（比如A的第i列和B的第j列的内积都算出来了），而循环里只算对应列的内积。浮点运算的顺序不一样，中间累积的舍入误差就会有差别——尤其是带权重的A'*W*B，运算步骤更多，误差更明显。如果要和循环结果完全对得上，得换个不做多余计算的向量化方法。

正确的向量化写法 #

无权重的情况（对应X1） #

用共轭转置的元素级乘法加求和，和循环的计算逻辑完全一致：

X3 = sum(conj(A) .* B, 1)';

这里conj(A)就是对A的每一列做共轭转置（对应循环里的A(:,ii)'），逐元素乘B之后按行求和，再转成列向量，结果和X1分毫不差。

带权重的情况（对应Y1） #

先算W*B，再和A做共轭元素级乘法后求和：

WB = W * B;
Y3 = sum(conj(A) .* WB, 1)';

或者直接写一行：

Y3 = sum(conj(A) .* (W*B), 1)';

这个步骤和循环里A(:,ii)'*W*B(:,ii)的运算顺序完全相同，数值结果和Y1没有差异。

最快的实现方式 #

几种方法的速度对比 #

• 循环：MATLAB的JIT现在对简单循环优化不错，但列数多（比如200）、行数大（2000+）时，循环的开销还是比向量化大。
• 矩阵乘法取对角：平白算一堆没用的非对角元素，纯纯浪费计算资源，速度肯定不如只算需要的部分。
• 元素级乘法+求和：只计算对应列的内积，没有多余运算，而且MATLAB对元素级运算和sum函数的优化非常到位，这是最快的方案。

额外提速技巧 #

如果你的权重矩阵W是Hermitian矩阵（也就是W == W'），可以保持写法不变，MATLAB会自动利用对称性优化。要是矩阵特别大，还可以考虑用GPU加速（如果有支持CUDA的显卡），但常规CPU环境下，上面的元素级乘法+求和就是最优解。

内容的提问来源于stack exchange，提问作者Robert Seifert