在 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，得到重复计算之一的值。

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

这是结果代码：

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秒，修改后的代码仅耗时0.4秒。这不是矢量化的好处，它只是预分配的好处，而不是一遍又一遍地重复相同的计算。

接下来，我们可以通过向量化运算来移除另外两个循环。这里的难点是 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 倍原代码。

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

Optimizing matrix calculations in for loops in Octave

loops

matrix

vectorization

octave