优化嵌套 for 循环以计算矩阵行的 xcorr

Question

我有 2 个嵌套循环，它们执行以下操作：

• 获取矩阵的两行
• 检查索引是否满足条件
• 如果这样做：计算两行之间的 xcorr 并将其放入新向量中
• 找到子向量最大值的索引，并将LAG矩阵的元素替换为该值

我不知道如何通过矢量化或其他方式加速此代码。

b=size(data,1);
F=size(data,2);
LAG= zeros(b,b);  

for i=1:b
    for j=1:b
        if j>i
            x=data(i,:);
            y=data(j,:);
            d=xcorr(x,y);
            d=d(:,F:(2*F)-1); 
            [M,I] = max(d);
            LAG(i,j)=I-1;
            d=xcorr(y,x);
            d=d(:,F:(2*F)-1);
            [M,I] = max(d);
            LAG(j,i)=I-1;
        end
    end
end

Answer 1

首先，关于浮点精度的说明...

您在评论中提到您的数据包含整数 0、1 和 2。因此，您希望 cross-correlation 给出整数结果。然而，由于计算是在double-precision 中完成的，因此似乎引入了一些floating-point error。此错误可能导致结果比整数值稍大或稍小。

由于您的计算涉及查找maxima 的位置，因此如果存在重复的最大整数值并增加了精度误差，您可能会得到略有不同的结果。例如，假设您希望值 10 为最大值，并出现在向量 d 的索引 2 和 4 中。您可能会以一种方式计算 d 并得到 d(2) = 10 和 d(4) = 10.00000000000001，但会增加一些精度误差。因此，最大值将位于索引 4。如果您使用不同的方法来计算 d，您可能会得到 d(2) = 10 和 d(4) = 9.99999999999999，错误方向相反，导致最大值位于索引中2.

解决方案？ Round首先是你的互相关数据：

d = round(xcorr(x, y));

这将消除浮点错误并为您提供您期望的整数结果。

现在，开始实际的解决方案...

解决方案 1：非循环选项

您可以将一个矩阵传递给xcorr，它将对列的每个成对组合执行互相关。使用它，您可以像这样完全放弃循环：

d = round(xcorr(data.'));
[~, I] = max(d(F:(2*F)-1,:), [], 1);
LAG = reshape(I-1, b, b).';

解决方案 2：改进循环选项

对于上述解决方案，data 的大小是有限制的，因为它会产生可能超过 maximum array size available 的大中间变量和输出变量。在这种情况下，for 循环可能是不可避免的，但您可以改进上面的 for 循环解决方案。具体来说，您可以为一对 (x, y) 计算一次互相关，然后只需翻转这对 (y, x) 的结果：

% Loop over rows:
for row = 1:b

  % Loop over upper matrix triangle:
  for col = (row+1):b

    % Cross-correlation for upper triangle:
    d = round(xcorr(data(row, :), data(col, :)));
    [~, I] = max(d(:, F:(2*F)-1));
    LAG(row, col) = I-1;

    % Cross-correlation for lower triangle:
    d = fliplr(d);
    [~, I] = max(d(:, F:(2*F)-1));
    LAG(col, row) = I-1;

  end

end