在 Octave 的 for 循环中优化矩阵计算

Question

我将代码从 Matlab 导入到 Octave，某些函数的速度似乎下降了。 我研究了矢量化，但以我有限的知识无法提出解决方案。 我想问一下，有没有办法加快这个速度？

    n = 181;
    N = 250;
    for i=1:n    
        for j=1:n
            par=0;
            for k=1:N;
               par=par+log2(1+(10.^(matrix1(j,i,matrix2(j,i))./10)./(matrix3(j,i).*double1+double2)));
            end
            resultingMatrix(j,i)=2.^((1/N).*par)-1;
        end
    end

尺寸在哪里：

    matrix1 = 181x181x2,  
    matrix2 = 181x181 -->  containing values either 1 or 2 only,  
    matrix3 = 181x181,  
    double1, double2 = just doubles

Answer 1

这是我的测试代码，我通过制作一些随机矩阵完成了你的代码：

n = 181;
N = 250;
matrix1 = rand(n,n,2);
matrix2 = randi(2,n,n);
matrix3 = rand(n,n);
double1 = 1;
double2 = 1;
tic
for i=1:n
   for j=1:n
      par=0;
      for k=1:N
         par=par+log2(1+(10.^(matrix1(j,i,matrix2(j,i))./10)./(matrix3(j,i).*double1+double2)));
      end
      resultingMatrix(j,i)=2.^((1/N).*par)-1;
   end
end
toc

请注意，k 循环内的代码不使用k。这使得循环是多余的。我们可以很容易地删除它。循环执行相同的计算 250 次，将结果相加，然后除以 250，得出其中一次重复计算的值。

另一个重要的事情是preallocateresultingMatrix，以避免它随着每次循环迭代而增长。

这是生成的代码：

tic
resultingMatrix2 = zeros(n,n);
for i=1:n
   for j=1:n
      par=log2(1+(10.^(matrix1(j,i,matrix2(j,i))./10)./(matrix3(j,i).*double1+double2)));
      resultingMatrix2(j,i)=2.^par-1;
   end
end
toc
max(abs((resultingMatrix(:)-resultingMatrix2(:))./resultingMatrix(:)))

最后一行计算最大相对差。在我的 Octave 版本中是 9.9424e-15。它会因版本、系统等而有所不同。此误差是浮点舍入误差。请注意，原始代码将相同的值相加 250 次，然后除以 250，会产生比修改后的代码更大的舍入误差。例如，

x = pi;
t = 0;
for i = 1:N
   t = t + x;
end;
t = t / N;
t-x

给出-8.4377e-15，与我们在上面看到的类似的舍入误差。

原代码耗时81.5 s，修改后的代码仅需0.4 s。这不是矢量化的增益，它只是预分配的增益，不需要一遍又一遍地重复相同的计算。

接下来，我们可以通过向量化操作来移除其他两个循环。这里的难点是matrix1(j,i,matrix2(j,i))。我们可以使用(1:n*n).' + (matrix2(:)-1)*(n*n) 生成每个n*n 线性索引。这不是微不足道的，我建议您考虑一下这种计算是如何工作的。您需要知道线性索引计数，从左上角数组元素的 1 开始，首先向下，然后向右，然后沿着第 3 维。所以1:n*n 只是按顺序排列的二维数组中每个元素的线性索引。如果我们需要访问第 3 个维度的第 2 个元素，我们会在其中添加n*n。

我们现在有了代码

tic
index = reshape((1:n*n).' + (matrix2(:)-1)*(n*n), n, n);
par = log2(1+(10.^(matrix1(index)./10)./(matrix3.*double1+double2)));
resultingMatrix3 = 2.^par-1;
toc
max(abs((resultingMatrix(:)-resultingMatrix3(:))./resultingMatrix(:)))

此代码产生的结果与我之前的版本完全相同，运行时间仅为 0.013 秒，比非矢量化代码快 30 倍，比原始代码快 6000 倍。 p>